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Abstract. We have recently written a new code to simulate the long term evolution of spherical clusters of stars. 
It is based on the pioneering Monte Carlo scheme proposed by Henon in the 70's. Unlike other implementations of 
this numerical method which were successfully used to investigate the dynamics of globular clusters, our code has 
been devised in the specific goal to treat dense galactic nuclei. In a previous paper, we described the basic version 
of our code which includes 2-body relaxation as the only physical process. In the present work, we go on and 
include further physical ingredients that are mostly pertinent to galactic nuclei, namely the presence of a central 
(growing) black hole (BH) and collisions between (main sequence) stars. Stars that venture too close to the BH are 
destroyed by the tidal field. We took particular care of this process because of its importance, both as a channel 
to feed the BH and a way to produce accretion flares from otherwise quiescent galactic nuclei. Collisions between 
stars have often been proposed as another mechanism to drive stellar matter into the central BH. Furthermore, 
non disruptive collisions may create peculiar stellar populations which are of great observational interest in the 
case of the central cluster of our Galaxy. To get the best handle on the role of this process in galactic nuclei, we 
include it with unpreceded realism through the use of a set of more than 10 000 collision simulations carried out 
with a SPH (Smoothed Particle Hydrodynamics) code. Stellar evolution has also been introduced in a simple way, 
similar to what has been done in previous dynamical simulations of galactic nuclei. To ensure that this physics 
is correctly simulated, we realized a variety of tests whose results are reported here. This unique code, featuring 
most important physical processes, allows million particle simulations, spanning a Hubble time, in a few CPU 
days on standard personal computers and provides a wealth of data only rivalized by A'^-body simulations. 
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1. Introduction 

This paper is the second part of the description of the code 
we have developed in the past few years in order to inves- 
tigate the long-term dynamics of dense galactic nuclei. In 
a first paper (Freitag & Benz 2001, hereafter paper I), we 



presented the basic version of this Monte Carlo (MC) code 
which deals with 2-body relaxation. In this article, we add 
flesh to this kernel by incorporating physical effects that 
are of particular interest and relevance for galactic nuclei. 
The structure of the paper is as follows. In Sec 



1.1 



we motivate our interest in the dynamics of galactic nu- 
clei through a short review of the history of the study of 
this field. We then proceed to describe the new physics 
incorporated in the code. The principles of the basic ver- 
sion of our MC code are presented in paper I with which 
the reader is advised to get familiar; we don't repeat this 
information here. Our treatment of stellar collisions is 
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treated in Sec. ||, tidal disruptions in Sec. ^, while further, 
more minor, additions and improvements are described 
in Sec. ^. A variety of test simulations are reported and 
discussed in Sec. ^j. Finally we summarize this work and 
propose future developments in Sec. ||. An appendix is 
added to expose how we build initial conditions for use 
with our code. 



1.1. Astrophysical motivation 

Only very few reviews have been written about the dy- 
namics of galactic nuclei (Gerhard 1994 is the only recent 
reference known to us), hence, we introduce our work by 
a summary of the history of this complex field. 

The theoretical study of the dynamics of galactic nu- 
clei was initiated in the 60's, mostly to investigate whether 
stellar collisions in extremely dense clusters could power 
the, then recently discovered, quasars (QSOs). In these 
early speculations, the presence of a central massive black 
hole (MBH) wasn't assumed. QSO observations required 
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the concentration of a huge amount of matter in a smaU 
volume and it was thought that a stellar system could ex- 
pel its angular momentum and contract more easily than 
a purely gaseous configuration, hence reaching densities 
such that highly energetic collisions between stars should 
be commonplace (Gold et al. 1965; von Hoerner 196^ ). 



elusion that, provided a significant fraction of the emitted 
gas is accreted, it dominates the feeding of the BH in sys- 
tems of moderate stellar density while collisions are still 
the main player in denser nuclei and that the full range of 
AGN and QSO luminosities can be attained without hav- 
ing recourse to an external source of gas. More recently, 



While Woltjer (1964) proposed that collisions themselves 
would be a strong source of optical radiation and radio- 
emitting energetic particles, others pointed out that these 
disruptive events should lead to the formation of a massive 



compact gas object (Ulam & Walden 1964), maintained 
at QS O luminosity either by further collisions ( ^pitzer"fc 
Saslaw 1966 ; Spitzer fc Stone 1967 ) or m assive star for 
mation and supernovae (SN) explosions ( ganders 1970| ) 



Growtji of high-mass stars through repeated mergers in a 1992|- ICombes 2001 



Rauch (1999) has considered the relativistic dynamics of 
a compact stellar cluster dominated by a central MBH in 
an AGN and concluded that collisions, most of which are 
grazing, produce only little gas but may efficiently replen- 
ish the loss-cone for tidal disruptions. 

In the past decade, gas-dynamical processes have been 
increasingly favored over stellar dynamics as the main 
source of fueling of AGN ( Shlosman et al. 1990| ; 3hlosman 



low-velocity coUisional cluster was proposed by Colgate 
(196?!)] as another way of forming SN-powered QSOs. 
Unfortunately, in none of these early studies, was the stel- 
lar dynamics treated in a realistic way, most authors hav- 
ing recourse to some extension of the evaporative model 
of globular clusters (see, e.g., ^pitzer 1987 ). In particular, 
the process of gravothermal collapse was not known and 
the role of mass segregation not properly recognized. 

Nearly all further studies accounted for the presence 
of a central MBH, an object more and more widely ac- 
cepted as necessary to explain QSOs and others Active 
Galactic Nuclei (AGN), and likely to be present in at least 
some normal present-day nuclei, as a relic of past activ- 



that. 



ity (Lynden-Bell 1969). In a relaxed cluster were stars are 
destroyed in the vicinity of the BH, presumably by tidal 
forces ( Hills 19"75| ), their (quasi-)stationary density dis- 
t ribution was predicted to be a power-la w, cx _R~^/* 



and references therein). It is argued 
to achieve the highest QSO luminosities, the ini- 
tial stellar cluster has to be so dense that its formation 
is problematic and would, most likely, require to concen- 
trate a large amount of gas in the galactic center any- 
way. Furthermore, whether most of the gas emitted by 
stars -either in the course of their normal evolution or 
through collisions- finds its way to the MBH is uncer- 
tain (see Sec. |6.2| ). However, it may have been overlooked 
that the effective stellar relaxation rate, and, hence BH 
fueling through tidal disruptions or direct horizon cross- 
ings, may be highly enhanced by small departures from 
the assumption of a smooth spherical potential. Such de- 
partures may be the presence of o rbiting cores or nuclear 



BHs of smaller ac creted galaxies ( Polnarcv fc Rees 1994 
[Zhao et al. 2002|) , or triaxiality ( |Norman fc Silk 1983|) 



(Peebles 1972 



Bahcall fc Wolf 1976| , |l977D . The tidal dis- 



which may survive in the vicinit y of the BH even if it is 
destroyed at intermediate scales ( Poon fc Merritt 2002| ).F] 



ruption rate is dominated by stars t hat are brought onto 



Lightman fc Shapiro 1977 




Young et al. 1977 


; Cohn fc 


Kulsrui 197S, see Sec. 
cient to feed a QSO-cla 
is so high that collisions 
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). This rate is probably insuffi- 
vIBH unless the stellar density 
as production should dominate 


(Young et al. 1977; 


Young 1977; 


Shields & Wheeler 1978; 



Frank 1978 ). The link between the earlier BH-free stellar 
dynamical models for AGN a nd these studies of MBH- 
cluster systems was traced by Begclman fc Rees (1978 ) 



Even though purely stellar dynamical processes are 
probably only secondary in feeding QSO-class MBHs, they 
may be efficient enough to grow few million solar masses 
objects from BHs with a mass of a few hundreds Mq. 
Furthermore, questions regarding the interplay between 
the stellar nucleus and a central MBH are more press- 
ing than ever, as observational evidences for the pres- 
ence of MBHs in most, if not all, bright galaxies, in- 



who shpwed that most very dense stellar systems will nat- ct al 2000 
urally evolve to form large BHs. 



sive rate ( 


Ferrarese et al. 2001; [Vlerritt fc Ferrarese 2001; 


Kormendy fc Gebhardt 2001 


; Genzel et al. 2000 




Ghez 



and references in paper I). In particular, tidal 

-4 ,,^-1 



In the 80's, self-consistent simulations of the evolu- 
tion of galactic nuclei appeared, that confirmed these con- 
clusions ( [McMillan et al. 198l| ; [Duncan fc Shapiro 198^ ). 



models were based on Fokkcr-Planck and Monte 
odes first developed to study globular clusters. A 



These 
Carlo 

serious' shortcommg of these works was to assume that col- 
lisions were completely disruptive. (David et al. 1987a ^ 
and Murphy et al. (1991[ ) improved on this by implement- 
ing som e extension of the simple semi-analytical prescrip- 
tion of gpitzer fc Saslaw (1966 ) to account for partial 
disruptions but the introduction of collisions into Fokker- 
Planck codes had to be done in a quite unrealistic way (see 
Sec. 5^). Stellar evolution was also included with the con- 



disruptions at a rate of order 10^ yr~ seem unavoidable 
for BHs less massive than a few 10^ M©, with the likely 
consequence of bringing back to active life an otherwise 
quiescent galactic nucleus (Hills 1975; Lidskii & Ozerno: 



1979; Rees 1988; Phinney 1989; Bembayfc West 1993; 3yei 



Ulmer 1999; Magorrian & Tremaine 1999). Ironically, 



while tidal disruptions are deemed too rare to be the main 
contributor to the growth of MBHs, even in very dense nu- 
clei, they are predicted in present-day normal nuclei with 



^ Unfortunately, such possibilities, although pointing to the 
importance of stellar dynamical processes, could only be in- 
troduced approximately in our code which relies on spherical 
symmetry. 
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a rate which is embarrassingly high in regard to the low 
luminosity of these objects, a fact that has been used to 



impose constraints on gas accretion models (Sanders & 



van Oc 



stcrom 1984 



Menou & Quataert 2001). Some flar 



ing events in the UV or X-ray band from the center of 
active and non-active galaxies have been tentatively in- 
terpreted as the accretional aftermath of tidal disruptions 
( preiner et al. 200(]| ; [Komossa 200l[ |Rcnzini 200l| , and 
references therein). But further conclusions have to await 
more complete stellar dynamical simulations, like the ones 
we propose to carry out with our code, and a better under- 
standing of the post-disruption accretion process in order 
to predict its observational signature (wavelength, inten- 
sity, duration, etc.) Beside the accretion flares, another 
promising observational consequence is predicted: the pro- 
duction of h ot, very bright, stellar cores of tidally stripped 
giant stars (Pi Stefano et al. 2001). 



As i TOcallod above, colliaiono wore firot thought to play 



a majo r role in very donoo nuclei modclo, either by feeding 



1.2. General approach and limitations 

As is clear from this introduction and was already stressed 
in Paper I, the physics of galactic nuclei is a very intricate 
problem, with dozen of physical processes or aspects that 
can potentially play a role and interfere with each other. 
Any really general and realistic approach would have to 
face too many computational challenges and unknowns 
concerning the physics, initial and limit conditions to be 
feasible at the present date. Various numerical methods 
have different limitations and require different simplifying 
assumptions which delineate the class of models that can 
be treated. 

For instance, it is increasingly recognized that galaxy 
merging is a common process in the universe and that 
such ev ents have deep imprint on the structure of gala ctic 
nuclei (|Nakano fc Makino 1999| ; [Merritt fc Cruz 2()0l|) . Of 
particular interest is the for mation and evolution of binary 
BHs formed in the process dBegelman et al. 198(]|; [Gould &| 
Rix 2000|; [Hemsendorf et al. 200l|; [Milosavljevic fc Merrit 



the M BH or by directly producing the AGN luminooity. 
The latter class of mod els, now incorporating a centra l 
BH, has been revived by Torricelli-Ciamponi et al. (2000 ) 
but should be examined in the light of a more refined 
treatment of stellar collisions and stellar dynamics. Even 
if they are not frequent enough to have a strong impact 
on the dynamics or BH fueling, collisions may have inter- 
esting observational consequences, by producing peculiar 
stellar populations, like blue stragglers ( ^ills et al. 2001 



2001; Yu 2002). Self-consistent simulation of these highly 



and references therein, in the context of globular clusters 
or destroying giant stars (Genzel et al. 1996; Alexander 
1999; liailey & Davies 1999), for instance 



In addition to the now almost 'classical' questions con- 
cerning tidal disruptions and collisions, the stellar dynam- 
ics of galactic nuclei is key in other processes of high ob- 
servational importance. An important example is capture 
of compact stars on relativistic orbits around the MBH. 
Through relaxation or collisions, a compact star may get 
on a very elongated orbit with such a small pericenter 
distance that emission of gravitational waves will drive 
further orbital evolution until the star plunges through 
the horizon of the MBH (Hils & Bender 1995; 3igurdsson 



& Rees 1997; Freitag 2001; [vanov 2001). As these waves. 



if successfully detected and analyzed, would be a direct 



probe to the space-time geometry near MBHs (Thorne 
1998 |; [flughes 200l4 |b|), such relativistic MBH-star bina- 



ries will be prime-interest sources for the future space- 



borne laser interferometer LISA (Danzmann 2000). This 
question and other ones to be mentioned in Sec. 3.2 are 



beyond the scope of this paper and the relevant physics 
are not included in the code described here (see, however, 
Freitag 2001, for our first results concerning the capture 



of compact objects). Nonetheless, they strongly motivate 
the need for detailed numerical models of the stellar dy- 
namics in the center-most parts of galaxies. 



dynamical episodes in the life of galactic nuclei can only 
be done with Ai"-body codes in which the orbits of N par- 
ticles are explicitly integrated for many dynamical times. 
However, such direct iV-body integrations are extraordi- 
narily CPU-demanding and, when various physical pro- 
cesses interplay whose relative importance depends on N , 
their results can not be safely scaled to ^ 10^ to repre- 
sent a real nucleus. Hence, even wi th cutting-ed ge special 
purpose computers like GRAPE-6 (Makino 2001), A^-body 
simulations can not follow the evolution of a galactic nu- 
cleus over a Hubble time if relaxation is appreciable. 

The A'' barrier can only be broken through by trad- 
ing realism for efficiency. This is done mainly through 
three core assumptions: (1) Restricted geometry: we 
assume that the nucleus is of perfect spherical symme- 
try. (2) Dynamical equilibrium: at any given time, the 
system is a solution to the coUisionless Boltzmann equa- 
tion ( Binney fc Tremaine 198'^ ). (3) Diffusive 2-body 
relaxation: the departures from a smooth gravitational 
potential which is stationary on dynamical time scales, are 
treated as a large number of uncorrelated 2-body hyper- 
bolic encounters leading to very small defiection angles. 
This is the base of the standard Chandrasekhar theory of 
relaxation ( phandrasekhar 1960 ). 

To our knowledge, assumptions (2) and (3), which un- 
derlie the derivation of the Fokkcr-Planck equation from 
the Boltzmann equation (Binney & Tremaine 1987), are 
shared by all methods aimed at simulating the relaxational 
evolution of stellar clusters and all of them also rely on 
spherical symmetry, with the exception of the code devel- 
oped by pinsel fc Spurzem (199"9| ) and |Kim et al. (2002D 
which allows overall cluster rotation (see paper I for a 
short review of these various methods). We have based 
our cod e on the Monte Car l o (MC ) scheme invented by 



Henon ( |l971b| ; |l971a| ; |l973| ; |l975|) . The reason for this 



choice, presented in detail in paper I, is basically that this 
algorithm offers the best balance between computational 
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efficiency, with CPU time scaling like A^pln(ciVp) where 
Np is the number of particles and c some constant, and 
the ease and realism with which physics beyond relax- 
ation, in particular stellar collisions, can be incorporated. 
Other codes stemming from Hcnon's scheme have been 
developed and very successfully adapted to the dynam- 



Og; p^regeau et al. 2002p but we are not aware of tion during these epi sodes ( ^ffmann fc HachndF2000 



ics of globular clusters ( Stodolkiewicz 1982, 1986; Giersz 
1998 |, fioOlj ; jJoshi et al. 2000|; [Walters et a!l. 20"00| ; |Joshi 
et al 

any previously published adaptation of this method to the 
realm of galactic nuclei. 

In principle, there is no other restriction on the im 
tial conditions for the cluster than conditions (1) and (2) 



ing mass segregation and/or anisotropy, makes it difficult 
to test the implications of these implicit assumptions. 

The evolution of galactic nuclei is thought to go 
through highly dynamical phases, most noticeably merg- 
ers with other nuclei predicted by popular hierarchical 
structure formation scenarios. It is often assumed that 
the central BHs formed as intermediate mass objects 
(100 — 1000 Mq) and grew mainly by luminous gas accre- 



Volonteri et al. 2002, and references therein) but the op- 



posite view, i.e. that MBHs formed at high redshifts in the 
core of only a small fraction of proto-galaxies and grew 
mostly by merging together, cannot be ruled out (Mcnoi; 



In practice, however, the code we use to build the initial 
cluster (see Appendix) is limited to systems for which the 
distribution function (DF) depends on the energy only 
and doesn't account for the presence of a BH at the cen- 
ter. The ffist restriction implies that we cannot consider 
systems that present initial velocity anisotropy or mass 
segregation. The second forces us to start with 'seed' cen- 
tral BHs, i.e. the BH has to be initially so light that its 
addition at the center of the nucleus doesn't noticeably 
perturb the dynamical equilibrium. These limitations cor- 
respond to the class of models that have been investigated 
in most previous studies. For instance, with the exception 
of Ranch (1999 ), all the self-consistent simulations cited 
above considered a seed BH which grows through accre- 
tion of stellar matter. Even though this is not a favored 
BH growth scenario anymore, in this paper, we adopt such 
models mainly as a mean to establish the correct working 
of our code through comparisons with the literature. 

If the central BH forms on a time scale much shorter 
than relaxation time but longer than dynamical time, pre- 
sumably through infall of gas from outside the nucleus, as 



et al. 2001). Anyway, although the MC code cannot follow 



proposed by, e.g., van der Marel (1999) and MacMillan & 



Henriksen (2002), the stellar cluster reacts adiabatically, 
a process our code can cope with, as demonstrated in 



Sec. 5.1 . This allows to create models at dynamical equilib- 



rium which contains a central BH of significant mass. Our 
procedure for creating initial conditions can be adapted to 
clusters with central BH for whic h the energy-dependen t 
DF is known, such as 7- models ( Tremaine et al. 1994 ). 
In recent simulations to be reported in further papers, we 
use these models to investigate the dynamics of present- 
day galactic nuclei. The aim of this approach is to gain 
information about the rate and characteristics of interest- 
ing events (collisions, tidal disruptions, captures ... ) in 
z ~ galaxies without trying to guess which are the high-z 
"initial" conditions. However, it is observationally, as well 
as theoretically, doubtful that nearby galactic nuclei are 
devoid of anisotropy or mass segregation^. However, the 
lack of published generalizations of the ry-models, includ- 



^ For instance, we find for a model of the central cluster of the 
Milky Way, that significant segregation of stellar BHs appears 
in less than 1 Gyr so that assuming that no segregation has 
occurred in the past history of the system is unrealistic. 



these dynamical phases, one can easily use the outcome 
of iV-body simulations of such processes as initial condi- 
tions, as soon as dynamical equilibrium is reached and the 
system is reasonably spherically symmetric. The explicit 
knowledge of the distribution function is not required here 
because one can directly turn each A^-body particle into 
one or a few super-star (s). 



1.3. Units and definitions 

When we do not explicitly indicate astrophysical units, we 
use the "code" units defined in Sec. 3.2 of paper I. 

G is the gravitational constant. A/bh is the mass of the 
central BH, Mc\ the total stellar mass and i?ci the radius 
of the cluster (if finite) . We use the following definition for 
the core radius: Rc = \J / AttG pq where ctq is the central 
ID velocity dispersion and po is the central density of the 
cluster. 

We assume the following relation for the Coulomb log- 
arithm: ln(7A'^) with 7 = 0.4 for single-mass models and 
7 = 0.01 when there is an extended stellar mass spectrum. 
A^^, is the total number of stars in the cluster. In prin- 
ciple, the Coulomb logarithm, A should be proportional 
to A'^, only if the cluster is self-gravitating. In a central 
region of radius GMbhct"^ ~ Rc\{MB}i/Mc\) (assuming 
A/bh ^ Afci, cTv is the velocity dispersion of the stars 
far from the BH) , the BH gravitationally out- weights the 
stellar cluster. There, the velocity dispersion at distance 
R of the center is cr1{R) — GM-qh/R and a steep cusp 
of stars is expected to develops so that, 6max — i? is a 
sensible choice. Consequently, according to Eq. 6 of pa- 
per I, A oc A/bh/AT* seems more appropriate (Bahcall & 
Wolf 19761; iLightman & Shapiro 19771; pVIiralda-Escude & 



Goidd 200c ). We have conducted test calculations with 
a _R- variable Coulomb ratio set to A cx rorb(^)/7min(-R) 
where Torb ~ {GMr/ R"^)~^^^ is a measure of the orbital 
time and Tmin corresponds to the shortest effective 2-body 
encounter, i.e. Tmin ~ bo /ay ~ GM^a~^. Such a choice is 
motivated by the fact that a transient potential fluctuation 
with time scale much longer than Torb will act adiabati- 
cally on the motion of a given star and thus leave its orbit 
unchanged after it is over. Results are not significantly af- 
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fected by the choice of A, which convinced us to keep the 
simple A = jN^ relation. 

2. Stellar collisions 

2.1. Use of SPH collision simulations 

The inclusion of realistic collision^ is probably the main 
improvement over previous cluster evolution codes that 
our scheme features. In the past few years, we have been 
computing thousands of 3D hydrodynamics simulations of 
collisions between MS stars using a SPH code dBenz 1990| ). 
For simplicity, only collisions between main sequence stars 
are considered for the time being. Actually, giant stars are 



expected to dominate the collision rate (Bailey & Davies 



1999; Freitag & Benz 2002a). The effects of collisions are 



included in the cluster simulations with unpreceded re- 
alism by interpolating the outcome of these events from 



the huge SPH-generated results database (Freitag & Benz 



2002a). In so doing, we get rid of many of the uncertain- 
ties introduced by the simplistic recipes formerly used in 
simulations of coUisional cluster dynamics. Unfortunately, 
even with such a procedure, important simplifications of 
the physics have still to be done. The major ones are con- 
nected with the possible formation of binaries through 
tidal dissipation of orbital energy and to the stellar evo- 
lution of the stars that have suffered from a collision. We 
discuss both problems in turn. 

The cross sections for the formation of so-called "tidal- 



binaries" are not well known (Press & Teukolsky 1977; 



Lee & Ostriker 1986; 


McMillan et al. 1987; 


Benz & Hilk 


1992; Lai et al. 1993| 


Kim & Lee 1999) and their long- 



term evolution is still debated (Benz & Hills 1992; Lai 
|et al. 19^ ). Hence, it is fortunate that the rate of tidal 
captures is overtaken by the rate of collisions as soon as 
(Ty/v^ > 0.1 where a-a is the ID velocity dispersion of the 
stars and = •\/2GM,/_R» is the escape velocity at the 
surface of a star (Lai et al. 1993, Fig. 16). Consequently, 
as we expect quite high stellar velocities in the center of 
galactic nuclei (particularly near a super-massive BH) , we 
decided to neglect tidal capture in our code. 

A parameter of prime importance is the star's radius 
as it determines its collisional cross section and, hence, 
the probability of subsequent collisions that could lead, 
for instance, to the runaway build-up of more and more 
massive stars by multiple mergers. After a collision, as a 
large amount of energy has been injected into the stellar 
envelope, the star is much larger than a MS star with the 
same mass. However, on a Kelvin-Helmholtz time scale 
(2kh) the radius shrinks back to the MS value, as the 
stellar structure returns to thermal equilibrium. Here we 
assume Tcou 3> Tkh so we can neglect the short swollen 
phase and attribute a MS radius to the collision product. 

When stellar evolution is taken into account, it be- 
comes in principle necessary to know what amount of 



Here, by "collision", we mean a genuine hydrodynamical 
contact encounter between two stars, as opposed to mere 2- 
body gravitational deflections. 



collisional mixing occurs and how it affects the MS life- 
time Tms of the product. We can expect that, contrary 
to parabolic mergers where only little mixing takes place 
(Sills et al. 1997, for instance), high velocity collisions are 
able to rejuvenate the star by bringing fresh hydrogen- 
rich gas from the outer parts to the center. If two stars 
of unequal masses merge together, simulations show that 
the smaller one, whose material is of lower entropy, sinks 



to the center of the larger one (Lombard! et al. 2002). 
This appears as an efficient mechanism to bring fuel di- 
rectly to the core of the large star and delay hydrogen ex- 
haustion. Conversely, the higher mean molecular weight /i 
that results from spreading the central Helium (produced 
by H-burning on the MS) leads to a important decrease 
of Tm s compared to a star with a "normal" composi- 
tion ( IClaret fc Gimenez 199^ ). I ndeed, fr om homological 
relations, one finds: Tms oc /i^** (Kippcnhahn & Weigert 



1994 ) . On the other hand, the radius depends only weakly 
on /i (i? cx f/^-^ for the CNO-cycle) so we can safely ne- 
glect the effects on the collision cross section in our simu- 
lationsQ. 

2.2. Collision rate. 

Let's consider a close approach between two stars with 
masses and radii Mi, Ri and M2, R2, respectively. The 
relative velocity at infinity is v^ci and the impact parame- 
ter b. A collision occurs when the centers of the stars are 
closer to each other than d = ry(i?i -|- R2) {ij = 1 for gen- 
uine collision, rj < 1 for merging, 77 > 1 for tidal capture 
when Wrci is small enough). Until this collision distance 
is reached, we neglect the gravitational influence of other 
stars as well as any mutual tidal interaction. So the prob- 
lem reduces to a simple hyperbolic approach between two 
point masses. This gives us, the largest impact parameter 
leading to contact, 6max, and the cross section, 

2 ^ 



S. 



(12) 
coU 



7:bi 



Trf]^(Ri+ R2y' 



Wrel 



where 



,(12) 



'2G'(A/i +M2) 
i?l -|- i?2 



(1) 



(2) 



The second term is the bracket of equation [I] is the gravi- 
tational focusing which enhances the cross-section at low 
f 12) 

velocity (5*^011 oc i?i -I- R2). At high velocities iScoii tends 
to the geometrical value iTp{Ri + R2Y ■ So, the collision 
rate for a test-star "1" in a field of stars "2" having all the 
same mass, radius and relative velocity to "1" is simply 

(i;2) 

= n2Vr<,\S)}2^ (3) 



diV< 



coll 



di 



^ How the outcome of further collisions will be influenced by 
structural changes due to previous collisions has not yet been 
assessed. This can be of importance in the case of "run-away" 
mergers. 
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where 712 is the (local) number density of stars "2" . 

If we are interested in the overall collision rate in a 
star cluster, the next step to do is to introduce a velocity 
distribution. Before considering more general cases, let's 
assume that all stars in the cluster have the same mass 
Af, and radius i?, (so we can drop over-score "(i^)" jj^ 
and ^coii) and that their density is n*. The average local 
collision time Tcoii{R) is found by integrating the collision 
rate (equation ^) over all possible velocities of the two 
stars. 



1 



coll 



— / d^Vid^V2f{vi)f{V2)\\vi~V2\\Scol\- 

n* J 



(4) 



As shown in Binney fc Tremainc (1987 ), the result for a 
Maxwellian distribution is 



1 



TcoU 



M 



= 16v^ r]^Rln^av ( 1 
28.4 



477(72 



(5) 



0,01 =- 




0,001 ^ 



0,000 



For the Plummer model, the result is very similar to Eq. ^, 
with only the numerical constant replaced by 28.6. 

The total number of collisions per unit time in the 
cluster is given by the integration of l/Tcoii(i?) over the 
whole cluster: 



Fig. 1. Collision rate as a function of radius in a Plummer 
cluster with 60 = 0.725 and = 10^. The solid line is 
the theoretical rate based on Eq. |^. The dots are statistics 
from a MC simulation run with no cluster evolution. "N- 
body units" are used (see Sec. Oh. 



diVeoll 



dt 



2.3. Relative collision rates between stars of different 
masses. 



For a Plummer model of total mass M, star number TV* 
and scale radius i?p, the collision rate by unit radius reads: 



diVcoii 
dt dR 



with 



(R) 



54^/2^-L X 



(l + u 



Rp el 

2N-21/4 



i + eo(i 



,1/2 



(7) 



R 3 M 

r;^ p" = p^^^-a^M^ 



and 



60 



v: ^ 3 Rp 

ArjaUO) VN. R. 



(Safronov number). 



As a check of our code. Fig. |i| depicts this rate along with 
the statistics produced in a inventory run during which 
the cluster's structure as a Plummer model was frozen. 

Carrying out the radial integration, we finally get the 
total collision rate in the whole cluster: 



dN, 



coll 



dt 



1 



Vg;^^ (4.25 + 5.2060) (8) 



We now address the case of a cluster with a distribution of 
stellar masses. For simplicity, we consider a discrete mass 
spectrum with A'sp components: Af, G {Mi}f^l with (lo- 
cal) densities n^. So, using Eq. |^, the rate by unit volume 
for collisions between stars of classes i and j, with veloci- 
ties Vi and Vj is 

r,,d^v,d''v, = Mv,)fjiv,)\\v, - v,\\S^^i^^,d^v,d''vj (9) 

where fi, fj are the phase-space DFs which are assumed 
to comply with (spatial) spherical symmetry and isotropy. 
Their i?-dependence is implicit. If we further assume 
Maxwellian velocity distributions with 1-D velocity dis- 
persions ffi and (Tj, the distribution of the relative ve- 
locity t>roi = Vi ~ Vj is Maxwellian too, with dispersion 

erf + 0-^. We keep Wroi 



||tJrci|| as the only rele- 
vant velocity variable by integrating Eq. M over the others: 



rcl 

cx riirij (Ri 



R. 



^2„,3 



t'rol 



rcl ■ 



(10) 



For a continuous mass spec- 
mass function as ip{Mi) = 



rj has been set to 1. 
trum, we define the 

n~^dn^{Mi)/d{logiQ{Mi)) so we have to substitute dru = 
n^^{Mi)d{logiQ{Mi)) for in the previous formula. In or- 
der to get a equation for the relative collision rate between 
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stars of different masses (per unit volume, log]^o(Mi) and 
\ogiQ{Mj)), we assume = aj = (T^V«,j and integrate 
over Vici- 

(11) 



with e(*^') 



4a2 " 2a2(i?, + i?,)- 



G{M, + Mj) 



For a Plummer model with no mass-segregation (and, 
thus, a unique (T„(i?)), this relation, when integrated over 
the whole cluster, leads to 

diVcoU 



rtot(M„M,0 



dtd(logioAf,)d(logioMj) 



Ri + i?i 



X ( 1 + 3.66 



, + Mj Rp 
~M i?i + Rj 

2 



(12) 



{R, + R,)/Rq ' ^ ^ 



with e = 3.66 



Rp/Rc- 

M/M(, 



In relation only the dependencies on stellar quantities 
have been preserved to insist on the relative collision rates 
between different stellar species. Although it relies on the 
admittedly unrealistic hypothesis of no mass segregation, 
Eq. |l^ proves useful as a prediction our code can easily 
(and successfully) be tested against, see Fig. ||. 

2.4. Introduction of stellar collisions in the MC code 

The difficulty of introducing stellar collisions in any stellar 
dynamics code is twofold. First, as the previous discussion 
has shown, it is not at all straightforward to determine 
the correct distribution of collision parameters (vrei, star 
types, position in the cluster, ... ). Secondly, provided 
the result of a particular collision is known (by perform- 
ing hydrodynamical simulations, for instance), we want to 
be able to preserve as much as possible of that valuable 
information when introducing it back in the cluster evolu- 
tion code. Due to their very structure^, some widely used 
schemes, based on an explicit resolution of the Fokker- 
Planck equation, impose such a highly simplified treat- 
ment of the collisions' outcome that it would not make 
much sense to devote energy to a realistic computation of 
these events. The MC method is exempt of such limita- 
tions. 



^ Their basic limitation lies in the principle they owe their ef- 
ficiency to: they model the stellar system as a set of continuous 
DFs (one for each different stellar mass). 



2.4.1. Global code modifications. 

Collisions introduce a new time scale in the code. There 
is consequently a new constraint on the time steps 



StiR) < fstfcoiiiR). 



(14) 



Tco\i{R) is an estimation of the local collision time. We 
chose the following definition, based on Eq. ||: 



Tcou(-R) 



= 16y/T:n^av{Ri 



(15) 



where cr^ = (^;2^ (bracketed) quantities are local aver- 
ages. This particular expression was chosen for its ease of 
evaluation and because, provided all stellar species have 
isothermal velocity distribution (quite a strong demand!), 
it reduces to exact relations in the two interesting limiting 
cases: 



coU 



C^coU/ 



S^Gn^-^^^ for « (vl) 
ay 

l6y^n^ay{Rl) for > (vl). 



(16) 

By requiring 

6m < fst {f-J{R) + fJ,{R)) . (17) 

we make sure that time steps are short enough to resolve 
both relaxational and collisional processes. Apart for this 
extended constraint, all the time-step determination and 
pair selecting machinery of paper I is left formally un- 
changed. 

2.4.2. Monte Carlo sampling of the collisions. 

Relaxation is due to the cumulative effects of a huge num- 
ber a small individual scatterings and can be treated as 
a continuous process, affecting progressively the particles' 
orbits. To be computationally tractable this phenomenon 
is discretized back into "super-encounters". In contrast, 
collisions do not act gradually but are genuinely discrete 
events, each of which strongly affect the properties of the 
implied stars. Hence, there seems to be no way to add up 
the effects of collisions into "super-collisions", no escape 
from the necessity to simulate them as individual events. 

When a pair of adjacent super-stars is selected to be 
evolved for a time step St, we randomly orient their veloc- 
ities and compute the local number density of stars of any 
kind, explained in paper I. The probability for a mu- 

tual collision to occur during that time span is, adapting 
Eq. I, 



,(12) 



(18) 



When compared to Eq. this expression could be thought 
to be an overestimate as n* is used instead of n2 ■ Actually, 
for a given super-star of type "1" , the expectation value 
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a) 



b) 




Fig. 2. Total relative collision rate 
rtot(-^^i, M2) between stars with 
masses Mi and M2 in a Plummer 
cluster without mass-segregation. 
The gravitational focusing param- 
eter is O — 1.5. Masses are in 
Mq . Lighter gray shades correspond 
to higher values. Successive contour 
levels correspond to factor of 2 de- 
crease in r. The same levels are 
drawn on both panels. The mass- 
function is ^'(M*) oc M~^-3^ for 
between 0.2 and 20Mq and the 
Mass-Radius relation is set accord- 



ing to stell ar m odels by pchallei 
et al. (1992| ) and Charbonncl ct al 



199£ ) a) Theoretical rate from 



Eq. |13|. b) Statistics from the MC 
code (400 000 collisions) with clus- 
ter evolution inhibited. This com- 
parison demonstrates the accuracy 
of the collision sampling in our code. 



for the number of collisions with super-stars of type "2" to sort out of this data and plug into the MC code. In an- 



is 



N, 



(12) 
coll 



■Prob. for neigh-" 
bor of being of 
.type 2 



Collision prob." 
if neighbor is of 
.type 2 



112 /n.j. 



P. 



(12) 



coll 



(19) 



as needed. The collision probability is compared with 

a random number Xiand with [0; 1 [-uniform deviate. If 

(12) 

-'^rand < ^coii ' ^ collisiou has to be simulated whose initial 
conditions are completely determined as soon as a value 
of the impact parameter b has been chosen according to 
probability density 



dP = 



2bdb ^ , , 

— if < £) < 6,nax, 

max 

otherwise. 



(20) 



Note that the super-star pair is tested for collisions before 
relaxation is applied to it. In case a collision is suffered, 
the orbits are probably deeply modified. So the relaxation 
step is skipped even if the pair survived. 

2.4.3. Treatment of an individual collision. 

As explained earlier, the outcome of collisions happen- 
ing in the course of the cluster's evolution is specified by 
a large set of 3D hydrodynamical simulations. These are 
potentially able to provide us with any detail, significant 
or not, about the state of the resulting star(s) and released 
gas. Most of this information, however is of no real rele- 
vance so we focus on the important parameters we have 



other paper (Freitag & Benz 2002a), we describe the way 
collisions are simulated with an SPH code and how we ex- 
tract the needed "macroscopic" information back from the 
simulation. Suffice to say that, if we assume the center of 
mass (CM) reference frames defined before and after the 
collision are the same (i.e. that M[w[ + Mi^w^ — where 
AI[ 2 and w'^ 2 are the post-collision masses and velocity 
vectors in the pre-collision CM frame), the kinematic out- 
come is entirely described by 4 numbers. They are M{, 
M2, the final relative velocity at infinity, 



2E' 



and the deflection angle ^coii- Further information is con- 
tained in the post-collision stellar structure but it may be 
ignored if one assumes, as we do, that the produced star(s) 
return to normal MS structure. These 4 numbers are all 
we need to implement collisions between super-stars fol- 
lowing exactly the same scheme as described in Sec. 4.2.1 
of paper I (steps 2-4) for purely gravitational encounters. 
The only added difficulty is connected with mass changes 
and the proper tracking of energy variation they imply. 

Note that when a collision between two super-stars oc- 
curs, it amounts to each star in the first super-star collid- 
ing with a star from the second super-star. As the number 
of stars per super-star is the same by construction, one can 
apply the outcome of the collision (new mass and velocity) 
uniformly to all stars of the super-star, i.e. to the super- 
star as a whole. When the stellar collision results in to 
surviving stars, we have to modify the orbital and stellar 
properties of both super-stars; when there is only one star 
left (merger or destruction of the smaller star, see Freitag 



M. Freitag & W. Benz: Monte Carlo Cluster Simulations II. 



9 



& Benz 2002b), one superstar is removed and the other 
one is given the properties of the remaining star; if both 
stars are destroyed, both super-stars are removed from the 
simulation. 

3. Tidal disruptions 

3.1. Loss cone theory 

If a star ventures very close to the BH, it may be bro- 
ken apart by tidal forces. The condition for an element of 
mass to be stripped away from the surface of the star is 
that the instantaneous gravitational attraction on it (due 
to the BH and the star itself) be lower than the required 
centripetal acceleration. In the simplified case of a non- 
rotating spherical star on a Kcplcrian orbit, this condi- 
tion determines the following disruption radius: 



A/, 



BH 



(21) 



\27r 

Where is the average density of the stellar matter. 
This approximation assumes Mbh ^ M^. Note that 
this is really only the condition for the tidal stripping 
of the outer layers of gas because the stellar density in- 
creases toward the center of the star. A more realistic ap- 
proach should account for elliptical or parabolic orbits, 
tidally induced deformation and the genuine hydrodynam- 
ical nature of this violent phenomenon. Moreover, if deep 
encounter certainly result in complete star destruction, 
milder ones would be responsible of partial envelope strip- 



ping. Many studies have addressed these aspects (Carter 



fe Luminet 1983; 


Evans & Kochanek 1989|; Laguna et al. 


1993; 


Fulbright 1996|; 


Ayal et al. 200(\). Fulbright per- 



formed SPH simulations of parabolic encounters whose 
strength can be parameterized by 



R 



pen 



For polytropic star models with n — 3/2 and n = 3, he 
found that stripping of half the stellar mass occurs for 
/3h — 0.8 and /3h — 1.7, respectively. In the present version 
of our code, complete disruption is assumed for /3 > /3h 
while the star is left undamaged for more distant encoun- 
ters. This corresponds to Eq. |l] with the factor 2^^^ re- 
placed by P^^. 

The "loss orbits" are the set of stellar orbits with 
pericenter distance -Rpori smaller than i?disr- For a star 
at distance R to the center with velocity modulus v, the 
loss cone (LC) is the set of velocity directions that leads 
i?a < ^disr, either going to the BH our coming from it 
(see Fig. ||). The aperture angle of the loss-cone, ^lc, is 
given by the relation 

GMbh a i?dii 



sm (t/LC 



) = 2 



/ -Rdisi 

V vR 



V 

y 



1 



R 



I ^.(^) ^4^disr) 



In case of a co-rotating spherical star on a circular orbit, 




disr 



orbital trajectory 



Fig. 3. Diagram of the loss cone. 



where $*(i?) = <I>(i?) -I- GMbh/R is the cluster contri- 
bution to the gravitational potential. As, for reasonable 
parameters, i?disr is a tiny value, typical loss orbits are 
very elongated, so that R ^ i?disr and GMBn/Rdisr — 
vUMbh/M^)^/^ > v^. Hence Eq. 1^ simphfies to 



q2 



^ GMBuRdu 



(24) 



The loss cone is usually very small, as is demonstrated by 
an order-of-magnitude estimate of 6'lc at the BH's "influ- 
ence radius" (i?i — GMbh/ct^): 



( Ah V'"^ R 



(25) 



2 X 10" 



( Mbh \ ^ / i?h \ ' 



for i?, = IRq and M, = IMq. 



(26) 



Rh is the cluster's half-mass radius. 

If it wasn't for relaxation or other orbit modifying 
mechanisms (collisions for instance) these loss orbits, if 
initially populated, would be drained over a dynamical 
time and no further tidal disruption would be expected in 
the subsequent cluster evolution, unless some increase of 
i?disr occurs. This could happen for the whole cluster as 
a result of the BH accreting gas supplied to it by other 
sources like stellar w inds (Mbh /"), or, as investigated by 
^yer fc Ulmer (199E ) , for those stars that experience rapid 
swelling when they become red giants {R^, /"). 

The crux of determining the rate of tidal disruptions, 
however, is the role of relaxation. This process is capable 
of replenishing loss cone orbits while at the same time it 
can remove stars from such orbits thus preventing them 
from being disrupted. These effects have been tackled ei- 
ther using quite rigorous approaches ( Lightman fc Shapirc 
1977| ; ICohn fc Kulsrud 1978| ; [Magorrian fc Tremainc 1999|) 



one get 5 a factor 3 instead of 2 inside (• • • ) 



1/3 



mainly aimed at their inclusion into Fokkcr-Plan ck codes, 
or resorting to more approximate descriptions (Frank fc 
Rees 1976t ^yer fc Ulmer 1999| ; [Miralda-Escude fc Gould 
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2000). Here we only outline the problem by recalling a few 
simple facts. 

Eq. 1^ can be recast in a simple characterization of 
loss orbits: 



J2 < J? 



LC 



2GAfBHi?d 



(27) 



a condition independent of energy E (for stars not too 
tightly bound to the BH). Thus the flux of stars to/from 
disruption orbits is chiefly controlled by J- "diffusion" in 
the vicinity of the Jlc borderline. For a given star, let 
SJorh be the mean quadratic variation of the angular mo- 
mentum due to relaxation during a single orbit (defined 
as the trajectory segment from a passage to apocenter po- 
sition to the next one), 



6Jc 



(28) 



If SJorh ^ Jhc, stars can survive many orbits, scat- 
tered into and out of loss trajectories before being tidally 
disrupted. It follows that orbits with J < Jlc are not 
strongly depleted and this regime is referred to as full loss 
cone. If the velocity distribution is initially isotropic, this 
process doesn't modify that fact and the fraction of stars 
disrupted per orbital period is simply those of velocity 
directions pointing in the loss-cone: 



dTVfuu 
diV 



9lc 



orb 



(29) 



Conversely, in the empty loss cone limit, SJorh ^ Jlc, 
there is no way back from the loss orbits and the situation 
can be described as a genuine diffusion in J-space. At a 
given energy, the star density in J-space gradually goes 
to zero as Jlc is approached from above. This negative 
gradient controls the diffusive flux of stars to the lethal loss 
orbits. Analytical treatment of this regime is far beyond 
the scope of this paper so we refer the interested reader 
to the above-mentioned previous studies and turn to a 
description of our MC approach to the problem. 

3.2. Implementation of loss cone effects 

A reliable determination of the tidal disruption rate re- 
quires for the numerical simulation of the relaxation pro- 
cess a resolution SJ^um < JhC in the empty loss-cone 
regime and SJnum = SJorh in the full loss-cone regime. 
The latter case could be treated by use of Eq. 



29 



as a 

quick shortcut but the former constraint cannot be cir- 
cumvented as easily. Unfortunately, whereas simulation of 
"normal" relaxation imposes a value of the numerical de- 
viation angle per step, (J^stop sufficiently smaller than tt 
(<5^stcp — '^/2\/Jst — O.Itt, see Eqs. 7 and of 10 paper I), 
resolution of the (empty) loss cone region is not attained 
unless (50stop < 9lg ^ tt! Furthermore, a foolproof ap- 
proach, not relying on a clear-cut a priori distinction be- 
tween "full" and "empty" regimes, would necessitate to re- 
duce SOstep to the tiny "elementary" orbital ^^orb step with 
a corresponding Stgtcp = Porh — ln(7iV*)A^^^Trei, thou- 
sands of times smaller than the desired (5istcp — fstTrci^- 



Although Shapiro (1985) was able to attribute such tiny St 
only to those particles orbiting close to (or inside) the LC, 
hence preventing too drastic a code slowing down, such a 
feature doesn't fit in any straightforward way into Henon's 
scheme. To mention but one impediment, the need of de- 
vising time steps that depend only on the super-star's ra- 
dial rank would impose St ~ Poib for a large fraction of 
super-stars. 

The simple structure of our code ~ mainly consisting in 
successive 2-super-star interaction steps - having proved 
to be both easy to grasp conceptually and reliable when 
applied to relaxational and coUisional simulations, we in- 
troduced loss cone effects in a way that required the least 
modifications. 

Let's consider a single step. If the encounter was a col- 
lision, we only need to test whether each surviving super- 
star entered the LC through the interaction and to disrupt 
it in such a case. Indeed, collisions are not to be refined 
into more elementary processes. On the other hand, after 
a gravitational super-encounter has been computed, with 
deflection angle J^stop in the encounter reference frame, 
each surviving super-star is examined for tidal disruption 
in turn by simulating its random walk (RW) in J-space 
during ^^stcp- In MC spirit, we estimate typical "represen- 
tative" for the diffusion angle during a single orbit, S9orh 
by scaling down SOgtcp to orbital time. 



-'orb 



%tcp with riorb 



St 



step 



(30) 



Let w be the super-star's velocity vector in the encounter 
frame. We decompose the step Stgtcp into a random walk 
of the tip of to on a sphere with fixed w — \\w\\ radius, 
starting at its initial direction. A brute force implementa- 
tion would require up to riorb steps of angular size S0orh, 
each one followed by a test for entry into the LC (Eq. p7| ). 
The number of orbits per Ststep typically ranging from 10"^ 
to 10^, such a procedure turns out to be extremely inef- 
ficient, requiring a huge number of operations to detect 
only a few tidal disruptions, even if super-stars with ini- 
tial velocities pointing too far from the LC are filtered 
out|^ Fortunately, the burden can be lighten enormously 
through use of adaptive RW steps. Indeed, n individual 
steps of length S with random relative orientation are sta- 
tistically nearly equivalent to a single "meta" one of length 
A = ^s long as A is sufficiently smaller than the 

distance to the LC, to keep the risk of missing a disrup- 
tion during these n RW steps at very low level. Here is the 
outline of the random walk procedure: 



Actually, as 5Sstcp oc \/ Ststep oc yfjst, the number of super- 
stars to be tested for entry into the LC per (mean) 5fstop scales 
roughly as (J^stop oc /^t, with norb oc f&t steps in each random 
walk. As the number of Jtstep needed to simulate the cluster's 
evolution for a given physical duration is oc f^t , the total 
number of RW steps scales as oc fst and the code gets slower 
for larger time steps! 

* More precisely, for planar RW, the length of the surrogate 
"meta" -step should be chosen according to a Gaussian distri- 
variance. 
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1. Preparation. The orbital period is integrated us- 
ing Gauss- Chebychev quadrature and (56'orb is deduced 
from Eq. 

2. Initialization. The initial angular coordinates (0, 9) 
oi w = {w^ , w''^ , w^) are computed. We set a variable 
L2 to the total quadratic deflection angle to be covered 
during SUtcp, L2 *- 

3. LC test. If w*g = + w^'f + (w^m + ^'"f < 
ULC = J'LC I R, the super-star has entered the loss cone 
and is disrupted. Otherwise, we proceed to the next 
step of the procedure. We recall that vqm is the ve- 
locity vector of the pair's center of mass in the cluster 
reference frame. It is considered constant during the 
RW process. 

4. Completion test. If L2 < 0, the random walk is over. 
We break from the RW loop, the super-star left unaf- 
fected. 

5. RW step. A new (meta-)step is realized. First its am- 
plitude is set according to 

A = max (^S0oih, min (^Amax, Asafc, v^) ) 7 (31) 

where A^ax — O.Itt and Agafo = Csafo(w*8 - wlc)/w 
with Csafo — 0.2-0.5. This relation ensures that meta- 
steps get progressively smaller, down to the "real" in- 
dividual (50orb when the loss cone region is approached 
during lo-RW. Then the (meta-)step direction on the 
sphere is set by an random angle, /3, with uniform 
[0, 27r[ deviate (see Fig. ||). This determines a new ori- 
entation ((/i, 6) for w. The remaining quadratic path 
length is updated, L2 ^ L2 — . The loop is closed 
by branching back to point ^ 

To conclude this section, we highlight some shortcom- 
ings in our treatment of the LC. Our procedure amounts 
to examining whether tidal disruption occurs during the 
fine-grained diffusion process numerically represented by 
a single super-encounter. Thus, as long as "normal", non- 
LC relaxation is concerned, the super-encounter and the 
explicit RW are two statistically equivalent descriptions 
of the particle's evolution during Stgtcp- But only if the 
RW process leads into the LC, is the particle's J mod- 
ified as this is needed to determine the outcome of the 
tidal interaction. Its energy isn't modified accordingly be- 
cause energy conservation would be violated if some en- 
ergy change were applied to the super-star without being 
balanced by an opposite modification for the other super- 
star that took part to the super-encounter^ The main risk 
is the introduction of some bias in the ^^-distribution of 
stars that endured partial tidal disruption. Furthermore, if 
the super-star survived the RW, we give it back the post- 
super-encounter orbital quantities. Hence, there a possi- 
bility that it will be left lying in the LC with no regards 

® Conversely, non-conservation of angular momentum 
doesn't show up explicitly for the contribution of any super- 
star to the total J is always zero, by spherical symmetry! 
However, there is a risk that such "hidden" non-conservations 
of J may reflect in the distribution of ellipticities by introduc- 
ing some nonphysical feature in it. 




Fig. 4. Geometry of one random walk step on the velocity 
sphere in the encounter reference frame. A^ is the adaptive 
i*^ step, Pi a random angle, Wi the particle's velocity after 
step i—1 and lOi+i its velocity after step i. Velocities with 
tangential component pointing in the shaded disk corre- 
spond to disruption orbits, i.e. with v^^ < wlc = Jlc/R 
in the cluster reference frame. v)jy^ is the tangential com- 
ponent of the pair's center of mass velocity. 

to its empty/loss nature! This means that the DF as rep- 
resented by the code is probably not accurate in the LC 
region. 

A possible cure to these problems would be to elim- 
inate the super-encounter phase and to perform a sym- 
metric RW for both super-stars at the same time. 
Unfortunately, this is not so easy for they do not share 
a common rioib- Also, consistency would dictate to start 
the random walk with the orbital properties of the super- 
star (which determine riorb, for instance) before the super- 
encounter. However, to save computing time, the RW's 
initial conditions are set to the orbital state modified by 
the super-encounter, as this spares an extra computation 
of the peri- and apocenter distances which are needed both 
to compute Porb and to select a radial position on the new 
orbit. Quite unexpectedly, tests have demonstrated that 
this trick doesn't introduce any significant change in the 
cluster's evolution (most notably, the BH's growth rate)p°[ 

In our description, we neglected the fact that if the 
BH is massive enough, its Schwarzschild radius i?s = 
2GMbh/c^ can exceed i?disr for stars with a given struc- 
ture so that they will be swallowed by crossing the horizon 

^" To be fair, the gain in speed is also quite modest, as most 
of computing time is spent in the orbital position selection 
procedure. 
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for tidal disruption, as a function of the 
mass of of the MS star, M^^^^'' ~ 1.6 x 

/ . , \ -1/2 / „ \ 3/2 / „ s -3/2 / ^ N -3/2 

lO^^o(^) (t) (^) (ll) . 
We assumed -Rpiungc = -^s and /3h = 0.8 and used Af»-i?, 
relations from realistic models of MS stars JSchaller 



et al. : 992; Mcynet et al. 1994; Charbonnel et al. 1999; 



Chabrier fc Baraffe 200C 



without being disrupted, i'or a star with solar mass and 
radius, this will happen for Mbh > 1.6 x lO^M© while, 
for giants with M = IMq and R = lOOi?©, only an un- 
realistic BH with Mbh > 1-6 x lO^^M© would be large 
enough to prevent disruptions from happening. In Fig. ||, 
we plot, as a function of the mass of the MS star, the max- 
imum BH mass for which tidal disruption can occur. Note, 
however, that assuming i?piungc = may be an underes- 
timate. Indeed, a particle with negligible energy at infinity 
would be pulled into the BH on a no-return in-spiral orbit 
by relativistic effects if its specific angular momentum is 
lower than J,nin = 4GAfBHC^^, as the effective potential 
does not have high enough a centrifugal rise. This critical 
J value corresponds to a parabolic orbit with pericenter 

def 

sep aration -Rpiunge = c?min = 4i?s in Newtonian mechan- 



ics (Shapiro fc Teukolsky 1983, Sec. 12.3). This should 



be used as the effective radius of the direct plunge sphere 
provided tidal disruptions occurring on such relativistic 
in-spiral orbits do not lead to observable accretion events, 
i.e. most of the stellar gas stays on in-spiral trajectories, 
which seems unlikely given the huge spre ad in orbita l en- 
ergy of the post-disruption gas elements ( Rees 1988 ). An 
easy modification of the code allows to account for direct 
plunges but, for the sake of comparison with results from 
the literature, they were not treated in any simulations 
presented here. 



In the present version of the code, we assume that each 
time a star enter the disruption sphere, it is completely 
shredded to gas and that all this gas is immediately ac- 
creted onto the central BH. Treating the accretion process 
as being instantaneous is certainly a good approximation 
when the mean time between successive disruptions (of or- 
der 10* years in present-day galaxies) is much longer than 
the time scale of individual accretion events (a few months 
to a few years). When this is not the case, one may as- 
sume that the gas piles up in some circum-BH reservoir, 
waiting to be accreted at a later time when the disruption 
rate has decreased and/or the increased B H mass allows 



a sho rter accretion time (see models a la Murphy et al 



1991| in Sec. |5^ 



On the other hand, assuming complete accretion prob- 
ably leads to an overestimate of the tidal feeding rate be- 
cause, due to the huge spread in the energy of debris, 
only 50 % of the stellar gas is left bound to the BH just 
after a complete tidal disruption (Rees 1988; Fulbright 



1996, among others). Furthermore, when the leading ex- 
tremity of this bound gas stream comes back to pericen- 
ter, it collides with slower moving material and shocks 
to such a high thermal energy that of order half of th e 
bound gas may eventually get unbound ( Ayal et al. 2000 ). 
Consequently, in future works, we should assume that only 
a fraction Saccr = 25 — 50 % of the tidally produced gas is 
accreted, but, to be consistent with other cluster simu- 
lations from the literature, all results reported here were 
obtained with eaccr = 100%. 

Finally, the assumption of complete disrupti on is also 
an ov er-simplification, as hinted to by, e.g., Fulbright 
(1996|) who showed that the transition regime between no 
damage and full disruption spans /? ~ 1 — !■ 3 for n = 3 
polytropes. Real MS stars with masses > IMq, not to 
mention giants, are even more concentrated than n = 3 
polytropes so that there is an important range of pericen- 
ter distances for which envelope striping, rather than com- 
plete disruption would result. Other non-disruptive tidal 
effects like spin- up ( Alexander fc Kumar 2001 ) are also 
of observational interest for the center of our Galaxy and 
we plan to extend the abilities of our code in order to be 
able to keep track of such "tidally perturbed" stars that 
can amount to an appreciable fraction of the inner stellar 
population (Alexander & Livio 2001). 



4. Other additions and improvements 

4.1. Stellar evolution 

Stellar evolution (SE) is, in principle, an important ingre- 
dient to incorporate in nuclei simulations. For a typical 
IMF, of order 40% of the Zero- Age MS (ZAMS) mass is 
lost from the stars in the first 10^° years, so SE is poten- 
tially one of the dominating source of fuel for the BH. Also, 
how stars are affected by relaxation, collisions and tidal 
disruptions obviously depends on their masses and radii. 
For example, compact remnants resist disruptive events 
and, with the help of mass segregation, may come to dom- 
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inate the central regions. Whether or not larger and larger 
stars may be formed through successive mergers also de- 
pends crucially on the relative time scales of stellar evo- 
lution and collisions. 

For the time being, our treatment of SE is simple- 
minded and straightforward. We assume that a star is 
"born" on the ZAMS and keeps the same mass and ra- 
dius during its MS life which is of duration Tms- W e 
use the relation TMs(Af*) given by Bressan et al. (1993| ). 
When it leaves the MS, this star is immediately turned 
into a compact remnant, according to the following pre- 



scription ( Miralda-Escudc & Gould 2000). All progeni- 
tors with masses lower than 8 Mq become 0.6 Mq white 
dwarfs, those with masses 8-30 become IAMq neu- 
tron stars and those with larger masses become 7 Mq BHs. 
Part of the emitted gas is accreted on the central MBH 
and the remaining is ejected from the cluster. This sim- 
plistic relation between the ZAMS mass of a star and the 
final product of its evolution mainly reflects the lack of a 
strong set of observational constraints or theoretical pre- 
dictions in this domain. In any case, it is known that the 
ZAMS — > remnant relation strongly depends on metallic- 
ity, if only because stellar winds do ( Maeder 1992 ). All in 
all, it appears to us that these aspects of SE are probably 
a main source of uncertainties affecting the prediction of 
stellar dynamical mechanisms in which remnants take an 
important part. 

SE introduces a new time scale, namely Tms in the 
present implementation. To resolve it correctly, we im- 
pose the time step 5t{R) to be smaller than a fraction 
fst (typically 0.05) of the minimum of Tms as evaluated 
in each cell of the same radial mesh we use to estimate 
Trci{R) and Tcoii(i?). But, contrary to relaxation and col- 
lisions, in the absence of a strong initial mass segregation, 
there is no reason for this time-scale to increase with in- 
creasing R. Consequently, when SE proceeds faster than 
other processes, it imposes (nearly) the same, very short 
St to all super-stars and we loose the advantage of in- 
dependent St. In the simulations we have performed so 
far with SE included, we assumed a unique initial episode 
of star formation a i = so that, as soon as high mass 
stars have been turned into remnants, the slowing down 
due to stellar evolution ceases and the total CPU time 
is only increased by a factor of a few. A more fearsome 
performance decline will result if some form of continuous 
stellar formation is simulated or if the red giant phase has 
to be resolved as well. 



4.2. Particle doubling 

To maintain a high resolution in the late evolutionary 
stages of a highly coUisional, disruptive or evaporative 
cluster, we resort to particle doubling. When the num- 
ber of remaining super-stars has reached half the initial 
number, every super-star is split into two copies with the 
same orbital and stellar properties. In the first stage of 
the procedure, both copies are left at the same position 



R where their "parent" was. Then, we pick each super- 
star in turn, in random order, and place it at a random 
position on its orbit, in a way identical to what is done 
at the end of a normal evolutionary step. In that way, we 
minimize the risk of maintaining potentially harmful cor- 
relations between super-stars descending from a common 
ancestor. Of course, after particle doubling, the number of 
stars represented by each super-star has to be divided by 
2. Some cluster models (like the one set according to DS82 
model E, see below) go through several episodes of parti- 
cle doubling. Implementing proper book-keeping was the 
main difficulty with this new, otherwise straightforward, 
feature. 

4.3. Miscellaneous 

Various minor improvements have also been recently 
added to the code. For instance, in order to ensure that 
the orbital parameters {E and J) and positions of the 
super-stars are given time to adapt to the (supposedly 
adiabatic, see Sec. |5.1| ) modification of the potential, we 
force time steps St{R) to be smaller than some fraction 

/cvap of the evaporation time, Tevap '= Mci(dMci/dt)~^ 
where Mc\ is the stellar mass of the cluster, and smaller 
than some fraction /i„t of the "intern mass evolution" time 

Ti„t(i?) =^ Mint(i?)(dMci/di)~i where Mint(i?) is the to- 
tal mass interior of R. Typically, values around 0.01 are 
used for /ovap and fint- 

Also, in addition to the usual test we perform each time 
a particle has to be evolved, we periodically check for all 
the super-stars to be bound. This is an iterative procedure 
because if, during the first pass, we detect super-stars that 
are unbound, we remove them from the system and this 
may unbound other particles. 

5. Test simulations 

5.1. Adiabatic adaptation of the star cluster to the 
growth of a central black hole 

In some instances the central BH can grow significantly on 
a time scale Tbh much longer than the cluster's dynamical 
time but still much shorter than relaxation time. Such a 
hierarchy naturally occurs if a substantial amount of gas 
is flowing into the BH from outside the nucleus. Quick 
BH growth can also happen if mass lost my stars, either 
due to normal stellar evolution (in a young cluster), or to 
disruptive collisions (in a very dense cluster) , is efficiently 
accreted on the BH. 

As a consequence of the slow modification of the 
potential, the shape of stellar orbits evolve while con- 
serving adiabatic invariants, i.e. the angular momentum 



J and the rad ial action In (Young 1980; Binncy & 
Tremaine 1987). Correspondingly, the density profile of 



stars around the BH and their velocity distribution are 
modified. Characteristics of the resulting stellar profiles 
have been worked out for various initial clusters, either 
semi-analytically, using the conservation of the DF when 
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Fig. 6. Adiabatic growth of a central BH in various cluster models. Evolution of the stellar density. Jagged solid lines 
are results of our MC simulations with 10^ super-stars. Smooth dashed lines are theoretical predictions based on the 
conservation of angul ar momentum and radial action. They have been computed with a code provided by G. Quinlan 
( Quinlan et al. 1995 ). The dot-dashed line segment indicates the asymptotic cusp slope from Eq. It applies for 



A/bh < -^^ci- a) Plummer model, b) Isochrone model, c) 7-model with 7 
between the MC results and the theoretical predictions is excellent. 



0. d) Hernquist model. The agreement 
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velocity anisotropy. Solid lines are our results, dashed lines are theoretical predictions from the code of Quinlan et al 
(1995|)] For the sake of display clarity, snapshots selected here are different from those in Fig. ^. Our curves have been 
smoothed with a sliding averaging procedure. To cover a larger range in radius, the average is done over a smaller 
number of super-stars at small and large radii than at intermediate positions. Given the high level of noise in the MC 
data, the agreement with Quinlan's predictions is very satisfactory until Mbh grows past O.SA/ci- From this time, the 
tangential anisotropy in the outer parts of our models fails to increase with larger BH masses (see text) . 



expr essed as a function of adiabatic invariants ( Young isotropy is conserved (Goodman & Binney 1984; Quinlan 



1980|; lt.ee fc Goodman 1989|; [Cipollina fc Berlin 1994 et al. 1995| ) 



Cipollina 1995 ; [Quinlan et al. 1995) or by means of 



A^-bod y simulation s ( Sigurdsson et al. 1995 ; Leeuwin fc 
Athan |ss"oula 2000|) 



These studies show that a power-law cusp develops 
inside the influence sphere of the BH, of radius i?,, in 
which GM-BYi/R exceeds the original velocity dispersion 
of the stars. According to (Quinlan et al. 1995), if the 
initial stellar cluster is isotropic, presents a density cusp 
p oc R~^^ with 7i > and a DF diverging near E = (/)(0) 
like f{E) (X {E — 0(O))~", then the final density cusp has 
an exponent 



7f 



7i 



7i 



(32) 



This result only applies very close to the BH if its mass 
is larger than the mass of the initial stellar core; the cusp 
may be steeper at intermediate distances, i?tians < R < 
R, with i?trans = R^ / R, ([Lcc fc Goodman 1989|; |Cipollina 



fc Berlin 1994, see also Leeuwin fc Athanassoula 200Cl| ) 



Another key feature is the development of noticeable tan- 
gential anisotropy in the central regions. In models with 
analytic cores (i.e. with (p(0) — p{R)) (x R^ near the cen- 
ter), this anisotropy, although it is caused by the central 
BH, does not actually appear in the center itself where 



We have performed simulations of the adiabatic growth 
of a central BH in a variety of cases. In addition to the tra- 
ditional Plummer model, we adopted the same set of mod- 
els as Quinlan et al. (1995). These are the isochrone clus- 
ter (|Hcnon 1959|, |1960|; [Binney fc Tremaine 1987|), which 



has an analytic core, and three '7-models' (Dehnen 1993 
[Tremaine et al. 1994[ ) whose density profile is 



PiiR) 



3-7 



Meli?b 



47r R-y{R + Rb) 



4—7 



(33) 



where i?b is the break radius. The used 7 values are 0, 1 
( Hernquist 1990| ) and 2 ( Jaffe 198^ ). None of these models 
has an analytic core. Eq. ^ predicts 7f = 3/2, 3/2, 2, 7/3 
and 5/2 for Plummer, isochrone, and 7 = 0, 1, 2 models, 
respectively (Quinlan et al. 1995). 

To simulate the process of adiabatic BH growth, we 
switched off relaxation and all the other physical processes 
in the MC code. The algorithm reduces then to moving 
super-stars on their orbits again and again (see Sec. 5.2 
of paper I) while Mbh is slowly increased. The time step 
condition is /int = 0.002 (see Sec. O). This relatively 
small value is required to get a correct evolution of the 
anisotropy in the outer parts of the cluster. With larger 
time steps, the particles at large radii react too impulsively 
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to the BH's growth and their orbits tend not to develop 
enough tangential anisotropy or even to become radially 
dominated. Note, however, that this problem only occurs 
when the BH's mass is larger than half the mass of the 
stellar cluster and that the density profile appears to be 
unaffected by this even for /int — 0.01. 

In Fig. |, we compare our results with the output of 
the code written by Quinlan et al. (1995) and kindly pro- 
vided by van der Marel. This code makes explicit use of 
the conservation of adiabatic invariants to determine the 
structure of the BH-embedding cluster and we can re- 
gard its results as secure predictions. As can be seen on 
these diagrams, the MC code behaves very nicely in this 
regime. Given the numerical noise to be expected from 
such a method, the density profiles are deemed to be in 
perfect agreement for all models. In Fig. 0, the evolution of 
the anisotropy profile for the Plummer and the Hcrnquist 
model is plotted. This quantity, when determined from 
MC results, suffers from a much higher statistical noise, 
so that a stronger smoothing must be applied to get useful 
curves. Despite this noise, it is quite clear that our results 
match the predictions very well, except for the outer parts 
that lack some tangential anisotropy for large A/bh, as al- 
ready discussed 
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Fig. 8. Growth of the central BH for models with ini- 



tial conditions similar to models B and E of Duncan & 



Shapiro (1982 ). Our results, obtained with 256k super- 



5.2. Cluster models with tidal disruptions 

Another idealized regime to which many theoretical and 
numerical studies have been devoted is the case of a the 
relaxed single-mass spherical stellar cluster with a cen- 
tral BH. Collisions are neglected but stars entering the 
tidal disruption re gion are destroyed and their mass is 
added to the BH. [BahcaU fc Wolf (1976| ) demonstrated 
that the quasi-steady state solution of the Fokker-Planck 
equation for this situation corresponds to a central den- 
sity cusp with p cx i?^^/*. Although these authors used an 
one-dimensional approach with the energy E as the only 
variable, more accurate numerical integrations of the sta- 
tionary FP equation in (E, J) space, with a proper account 
of loss cone effects, have confirmed this result (Lightman 
fc Shapiro 1977|; |Cohn fc Kulsrud 1978| ), as did evolution- 
ary models ( Duncan fc Shapiro 1983| , hereafter DS83, for 
instance). As testified by Fig. |1C|, we reproduce this re- 
sult. This plot shows the evolution of a model with same 
physics and initial parameters as model I of DS83 and is 
described in more details in the next sub-section. 

A few evolutionary models have been published that 
are based on these simple physical assumptions. Most were 
meant to explore the possibility of forming a MBH during 
the core collapse of a globular clustcip]. They usually start 
with a seed black hole which grows by consuming stars. 



However, neglecting the role of a mass spectrum and binary 
stars, they fall short of physical realism. Unless the cluster is 
born with a very high velocity dispersion, oc M/R K)rb) 
where M is the total mass, R a measure of the size of the clus- 
ter and Voih a typical value for the (internal) orbital velocity 
of binaries (an unrealistic assumption for globular cluster but 
which may apply to models of proto-nuclei of galaxy like those 



stars, (solid and dash-dotted lines) are compared with 
those of these authors (dashed lines). We made two sim- 
ulations of model E. Both have been stopped when the 
stellar cluster was reduced to 500 M©. In the first one 
(solid line), we start abruptly with a tidal radius smaller 
than the cluster which rapidly adapt to this truncation. In 
the second run (dash-dotted line), we let the cluster adapt 
gently to the tidal truncation before we actually start the 
simulation by switching on relaxation (see text). 



To check the tidal disruption rate given by our code, we 
compare the growth of the central BH in such models with 
results from the literature. 

Fig. PI shows such a com parison for models B and E 
of Dimcan fc Shapiro (1982 , hereafter DS82) to whom 
we refer for the specification of initial and boundary con- 
ditions. We have used the same setting as these authors 
except that, in our computations, there is no initial stel- 
lar cusp around the BH and that, for model B, the BH 
is present from the beginning of the simulation and not 
added at a later time as done in DS82. We don't think 
these minor changes have any significant effect because 
the initial BH amounts to only a tiny fraction of the clus- 
ter's mass (Mbh = 150, 250 Mq, respectively, with M^i = 
3 X 10^ Mq). The match between our results and those of 
DS82 is not very good. In particular, for model B, the BH's 



of Quinlan & Shapiro (199C)), the binaries will delay collapse 



and probably trigger core rebound before the central density 
is high enough f or efficient "tidal feeding " of a seed BH (see 



Giersz 2001 



Rasic 



Gao et al 199l| ; [Giersz fc Spurzem 2000 
3t al. 2001 , for simulations of globular clusters with primordial 
binaries). 
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growth starts at a significantly later time but produce an 
object of comparable mass. However, DS82's simulations 
were stopped shortly after core rebound, which does not 
allow a comparison at late times. Note that the growth 
starts when core collapse is sufficiently deep to bring many 
stars close enough to the BH to be disrupted and that it is 
stopped by the fact that the disruption of these stars, most 
of which have large negative energies, amounts to heating 
the stellar cluster. Consequently, the temporal shift be- 
tween DS82's growth curve and ours mostly reflects that 
our code predicts a longer core-collapse time, Tec. We re- 
fer to paper I for a discussion of this point and the large 
spread found in the literature for the value of Tec. 

Concerning model E, on the one hand, our value for 
the time of strongest growth, again a quantity nearly co- 
incident with Tec, nicely agrees with DS82. Note that 
this cluster, being a Plummer with a strong tidal trun- 
cation, evolves quicker and differently than an isolated 
cluster, which gives more weight to this agreement. At 
the end of our simulations, around 20Gyrs, the cluster 
has nearly completely evaporated. On the other hand, the 
BH's growth is steeper and stronger in DS82's simula- 
tion. There is no doubt that it would have produced a. 



signifirlantly larger final RH than in our ca,se. had their 
simulation been carried on up to cluster dissolution. The 
reason for this disagreement is not known to us. We sus- 
pected that it may be linked to the fact that, in our sim- 
ulation, the remaining cluster mass is lower at all times 
than in DS82, which may, in turn, be due to the way our 
and DS82's code cope with the strongly out-of-equilibrium 
initial conditions. Indeed ~ 10% of all stars are initially 
beyond tidal radius. In our model, the cluster loses 17% 
of its super-stars very quickly to adjust to the tidal trun- 
cation. To have a better handle on this problem, we re- 
made the simulation with a cluster model which was first 
allowed to settle to equilibrium with its tidal truncation. 
To do this, we "evolved" it with no relaxation or any other 
physical process but still moving super-stars on their or- 
bits in the usual way. If a selected super-star was found 
with apocenter beyond tidal radius, it had only a small 
proba bility (around 0.01) to be removed at this step and 



was ot lierwise kept (at the same position). We think that 
this method produces a better initial structure in which 
each super-star has been given time to react "adiabat- 
ically" to the enforcement of the tidal truncation. 15% 
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Fig. 9. Growth of the central BH for models with initial 
condit ions identical to those of Amaro-Seoane fc Spurzen^ 



(200l|) . Our results, obtained with 256k super-stars, (solid 



lines) are compared with those of these authors (dashed 
lines) . 



from a cluster with a relatively low number of stars can 
not be regarded as instantaneous: it takes of order one 
orbital time for a star to actually leave the cluster and 
the probability for it to be back-scattered onto a bound 
orbit is non vanishing. Whether or not some improvement 
in the line of this in our evaporation prescription would 
lead to a better agreement with DS82 concerning MBH(i) 
is not obvious as these two aspects may well be uncou- 
pled. Note that a similar mismatch in the BH's growth 
curve appears in comparisons with preliminary si mula- 
tions realized by Amaro-Seoane fc Spurzem (2001 ) with 
a gas code (see below), but doesn't s how up in compar- 
isons with other results obtained by Duncan & Shapirc 
(19831) with their MC code and by [Murphy et al. (199l| ) 



with a direct Fokker-Planck scheme (see next subsection). 
In Fig. 0, we display the growth of the central BH 



for clusters corresponding to the models used by Amaro- 



of the cluster mass is lost in this procedure and the re- Seoane & Spurzem (2001 



suiting cluster also shows less evaporation during its fur- 
ther, relaxation-driven, evolution. However, this does only 
increase the discrepancy with DS82 concerning the final 
mass of the BH, see the dash-dotted curve on Fig. ||. 
Our higher evaporation rate is probably due to our sim- 
pler prescription for escape. We immediately remove any 
super-star which gets on an orbit with apocenter distance 
beyond tidal radius, regardless of its actual position on 
this orbit. More realistically, DS82 allowed stars on escape 
orbits to be kicked back to bound orbits. Recent works 



( [Fukushige fc Heggie 200C| ; p^akahashi fc Portegies Zwart 
2000; Paumgardt 2001 ) made it clear that evaporation 



hereafter ASOl). These con- 
sist of 10^ 1 Mq stars distributed according to a Plummer 
density law with a core radius of 0.707pc. The cluster is 
seeded by a fixed central BH with an initial mass of 5, 
50 or 500 M0. Only the last 2 values have been used by 
ASOl. It is clear that for masses as low as 5 or even 50 Mq, 
neglecting the motion of the BH is quite an unphysical as- 
sumption which is required by the present limitations of 
numerical codes. Even if a close agreement is not reached, 
our results are very similar to the curves from ASOl. In 
particular, we get the same phenomenon of convergence 
at late times toward an unique value of Mbh- This value 
is however smaller by a factor of ^ 2 than that of ASOl. 
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Fig. 10. Evolution of the density profile for a cluster with 
initial conditions identical to models I/II of DS83. The ini- 
tial number of super-stars is 10^. As in DS83's model I, 
collisions are not simulated. One notes the rapid develop- 
ment of a p c>c i?"^ *"^ cusp. 



5.3. Galactic nucleus models including collisions 

After having checked individual aspects of the MC code 
in simplified models (pure relaxation in paper I, colli- 
sions rates in sections and 2.3, adiabatic BH growth 
in Sec. 



5.1 



), we turned to the few published works 
addressing the long term evolution of dense galactic nu- 
clei in order to check our code's global behavior in physical 
regimes more relevant to our astrophysical field of interest. 

We first wanted to avoid the extra complication of stel- 
lar evolution and discarded those papers which take it 
into account. Furthermore, by their nature, Fokker-Planck 
methods can only include coUisional effects in an approx- 
imate way so that they don't allow a clear ch eck of this 
aspect of the code. Finally, A^-body simulations ( Arabadjis 
1997 ; Ranch 1999 ), although much more realisti^^l, were 
deemed too noisy to provide reliable data to co mpare with 
So we chose the venerable models by Duncan fc 



Shapiro (1983 



hereafter DS83) to conduct tests that in- 
clude relaxation, tidal disruptions and stellar collisions. 
DS83 studied three different models. The initial structure 
is a King cluster with Wq = 8 made of identical stars with 



ditions 



A/©. Models I and II sharC the skmd initial Con- 



3.6 X 10 stars and a COre radius Rc = 0.50 pC 
(the to'tal radius is 34.7 pC). A s66d black hole is present 
at the center with an initial mass AfBH(O) = 5 x 10^ Mq. 
Model III was devised to reach quasar-like accretion rates. 



In the case of Arabadjis (1997 ), it is not clear, however, how 
reliably relaxation processes can be simulated with a TREE 
algorithm. 
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Fig. 11. Evolution of the density profile for a cluster with 
initial conditions identical to models I/II of DS83. The ini- 
tial number of super-stars is 2 x 10^. As in DS83's model 
II, collisions are simulated. They are assumed to be com- 
pletely disruptive. Instead of a steep p cx power 
law, the cusp in the center gets milder and milder. 



It initially contains 57 x 10* stars, it has Rc = 0.82 pc and 
A/bh(0) = 2 X 10^ Af©. Models II and III include stellar 
collisions. They are assumed to be completely disruptive 
and the gas they release is instantaneously and completely 
accreted on the BH. We used the same initial conditions 
and physics but, to assess the influence of the assump- 
tion of complete coUisional destruction, we carried out two 
extra simulations using our realistic, SPH-generated, pre- 
scriptions (models lib and Illb). 

In figures |l^ and |l^, we present the evolution of the 
density profile for models I and II, respectively. The most 
conspicuous feature of the first figure is a spreading central 
cusp with p cx Such a power-law profile is repro- 

duced here for the first time by a Henon-like Monte Carlo 
method. Fig. |ll| shows that when disruptive collisions are 
introduced in our calculations, as in DS83's model II, a 
much milder cusp first appears (with exponent ~ — 1) and 
progressively flattens (with exponent > —0.5). It has been 
repeatedly reported that collisions strongly decrease the 



steepness of the inner density profile (Duncan & Shapirc 
1983|; [Murphy et al. 1991| ; pavid et al. 19874 |b|; [Ranch 



1999). A slope of ^ —0.5 is often obtained. However, the 



simulations by Ranch (1999) point to the establishment 
of a flat, cusp-less central region, not unlike our own re- 
sults. Murphy et al. (1991 ) get a strong depletion of stars 
in the innermost part of the cluster, a result which is ap- 
parently reproduced in some of Ranch's models. For lack 
of resolution, there is no similar effect to be seen in our 
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Fig. 12. Evolution of the growth rate of the central BH in 
clusters with initial conditions identical to models I, II and 
III of DS83. Dot-dashed hnes are from DS83. Model I does 
not include stellar collisions. Models II and III treat them 
as causing complete disruption of stars. Solid lines with 
dots are our results for these systems. Dashed lines with 
dots (labeled lib and Illb) show the effects of a realistic, 
SPH-generated, prescription for the outcome of collisions 
which allows partial disruptions and mergers (see text). 
We used 512 000 to 2 x 10® super-stars in our simulations. 
'W-body" units are used. For models I and II, the time 
unit is 1.37 x 10^^ yrs and the unit for dM/dt is 2.6 x 
10~^ M0yr~^. In model III, these units are 9.81 x 10^^ yrs 
and 5.8 x lO'^Mgyr-^ 



simulations. The practical relevance of this discrepancy is 
probably low, however, because the size of this rarefied 
zone is so small that it would contain only a few Mq in 
most cases even without depletion. So the validity of a 
statistical treatment of such a tiny region is highly ques- 
tionable anyway. The evolution of the density profile for 
model III is qualitatively similar. Interestingly, model lib, 
which incorporate realistic, partially disruptive collisions 
also forms a cusp, but in the much denser model 

Illb, collisions are efficient enough to reduce the exponent 
to a value between -1 and -0.5. 



The growth rate of the BH is depicted in Fig. 12, The 
qualitative agreement with DS83 is satisfying even though 
the rate we obtain is higher by a factor of ~ 2 in initial 
phases of collisional models. The reason of this difference 
is unknown to us. The most important effect of a realis- 
tic treatment of collisional outcome is a strongly reduced 
accretion rate. This is mainly due to the fact that most 
collisions are grazing and consequently produce low mass 
losses even for high relative velocities. Indeed, neglecting 
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Fig. 13. Cumulative distribution of the fractional mass 
losses in collisions for a simulation of model II with im- 
proved treatment of collisions (see text). All collisions oc- 
curring before time T = O.lUt — 1.37 x 10^° yrs are in- 
cluded in this count. The solid line shows the number frac- 
tion of all collisions which resulted in a fractional mass loss 
lower than a given amount r). The dashed line indicates 
what mass fraction of coUisionally released gas came from 
collisions with fractional mass loss lower than rj. 



gravitational focusing, we get 



dN, 



coll 



ddn 



R1+R2 



for dniin < Rl + R2 



where dmin is the closest encounter distance for the equiva- 
lent 2 point-mass problem. The cumulative distribution of 
the fractional mass loss for model II is depicted in Fig. ^ 
Actually, the average mass loss per collision is as low as 
0.08 Mq despite an average relative velocity for collisions 
of Vrci = 8.8 (see Eq. ||). These examples clearly demon- 
strate that any incorporation of collisions in galactic nuclei 
dynamics must account for partially disruptive events. 

To conclude this series of tests, we turn to one of the 
most complete and widely used set of simulations of the 
long-term evolution of dense galactic nuclei published to 
date, namely the "d irect" Fokker-Planck integrations by 
Murphy et al. (199l| , hereafter MCD91). These authors 



included the following physics in their computations: 

— 2-body relaxation. It is treated in the standard 
Fokker-Planck way ( for a description of the m ulti-mass 
FP scheme see, e.g. Murphy fc Cohn 1988 and refer- 



ences therein). Note that, in the FP scheme, the cluster 
is represented as a set of DFs, each of which represent 
a discretized mass class, i.e., stars that have all the 
same stellar mass. 
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Stellar collisions. To get the mass loss for individual 
collisions, MCD91 use a semi- analytical method de- 
rived from the procedure invented by ^pitzer fc Saslaw 
(19^^. It works by decomposing the stars into thin 
columns of gas parallel to the relative velocity and im- 
posing conservation of momentum for each, completely 
inelastic, collision between a column from one star and 
the corresponding column of the other star. No lateral 
mass, energy or momentum transport is considered. 
The MS stars are assumed to be n = 3 polytropes 
with cx i?*. These mass- loss rates are then aver- 
aged over impact parameter and relative velocities to 
get rates that depend only on velocity dispersion and 
mass ratio which allows the authors to compute the 
instantaneous mass-loss rate for any mass class, due 
to collisions with stars from any other (or same) mass 
class. The total mass loss for a given time step and 
mass class is then converted into a number of stars 
to be removed from the class. This is obviously quite 
an inaccurate representation of the real way collisions 
change the masses of individual stars. Mergers are not 
included in this formalism. 

Tidal disruptions. Stars that get closer to the BH 
than the tidal disruption radius are assumed to be 
completely disrupted and their mass is instantaneously 
and fully accreted by the BH. Although our numerical 
scheme is widely different, we use basically the same as 



sumptions, here. Hence, we refer to MCD91 and Cohn 



& Kulsrud (1978) for a description of how this is im- 
plemented in FP codes. 
— Stellar evolution. A simple prescription is used in 
which stars stay on the MS for Tms(A/*) and then turn 
abruptly into compact remnants (CR). No giant phase 
is simulated and all mass loss occurs at the end of the 
MS. See MCD91 for the specification of Tms(^4) and 
the MS CR relation. 

The initial stellar clusters are Plummer models with 
a core radius of 1 pc. The total stellar mass is initially 
8.291 X (10^,108, 10^, 10*^)7^0 for models of classes "1", 
"2" , "3" and "4" , respectively. The stars are initially on 
the MS and obey a power-law mass spectrum, dA^^/dM* cx 
between 0.3 and 30 M©, with a = 1.5, 2.5 and 3.5 
for cases "A", "B", "C". 

The cluster is seeded with a BH of mass 
Mbh = IO'^Mq at its center. The BH eventu- 
ally swallows all the gas lost by stars, through 
normal evolution, collisions or tidal disruptions, 
but its growth rate is limited by the Eddington 
rate Me = LE/{r]c^) = 4:nGfieMBKmp/{r]caT) — 
2.5 X lO-^Moyr-i (77/0.1)-1(Mbh/10^Mo) where r/ is 
the efficiency factor for conversion of mass into radia- 
tion during the accretion process, /Xg is the molecular 
weight per free electron of the accreted gas (~ 1.13 for 
solar composition), nip the mass of the proton and ctt 
Thomson's cross-section. A "standard" value of = O.l is 
used. If the instantaneous rate of gas production from the 
stars, Mprod, exceeds Me, only an amount Me accretes 
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Fig. 14. Final BH mass for all the MCD91-like models. 
The lines connect models with the same IMF slope. We 
compare our results (dashed lines) to those from MCD91 
(dotted lines). Solid dots are for simulations with 256 000 
super-stars; the open star symbols are for B models with 
10^ super-stars. The triangles on the left axis indicate the 
total fractional mass loss due to stellar evolution for IMF 
with a — 1.5,2.5,3.5, at an age of 15Gyrs. This corre- 
sponds to the final BH's mass expected if stellar evolution 
was the only feeding process and no star could escape the 
nucleus. Letters indicate the process whose contribution 
to the final BH's mass dominates: "E" stands for stellar 
evolution, "C" for collisions and "D" for tidal disruptions. 
Most of the discrepancies between our results and those of 
MCD91 is due to the lower contribution of collisions (see 
text). 



on the BH while the remaining accumulates into a central 
"reservoir" -presumably an accretion disk- to be accreted 
later when Mpi-od has declined below Me. The gas is 
assumed to be funneled completely and instantaneously 
to the center, i.e. no gas remains in the stellar cluster 
or is expelled from the nucleus. The structure of this 
reservoir is not resolved in the simulations. Instead, it is 
assumed to be small enough to contribute to the potential 
as a central point mass, exactly as the BH. However, 
distributing the central mass in two components, the 
BH and this reservoir, can still influence the dynamics 
slightly through the fact that only the mass of the BH is 
used to compute the tidal disruption radius. On the other 
hand, interactions between the gas reservoir and stars are 
neglected (see Sec. |6.2| ). 

We have simulated all models specified by MCD91 with 
256 000 super-stars. For models of class B, we have redone 
the simulations with 10^ super-stars. We basically mimic 
the initial conditions and physics of MCD91. For instance. 
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Fig. 15. Evolution of the central mass (BH + gas reservoir) for MCD91-like models of class A (a = 1.5 for the IMF). 
The various hatching styles indicate the origin of the gas. The initial BH mass is too small to be visible on these 
diagrams (dark gray hatching). The thick line is the mass of the central BH, as limited by Eddington luminosity. Our 
simulations were realized with 256 000 super-stars. Note that the ordinate mass units are different in each panel. For 
this top-heavy stellar spectrum, the role of stellar evolution is clearly dominant even in model A where the high stellar 
density boosts the collision rate. Panels (a) to (d) correspond to decreasing initial cluster mass (see text). 
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Fig. 16. Same as Fig. ^ but for models of class B (a = 2.5). 
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Fig. 17. Same as Fig. |l^ and |l^, but for models of class C (a = 3.5). In this model with a stellar IMF strongly 
dominated by low masses, the role of stellar evolution is minimized so that collisions and tidal disruptions dominate 
the gas production rate. 
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Fig. 18. Evolution of the gas production rate for galactic nucleus models with initial conditions corresponding to 
models 1B-4B of MCD91. We plot the amount of gas the stars release per year through different channels: stellar 
evolution, collisions and tidal disruptions. Note that, at early times, only a fraction of this gas is accreted by the central 
BH while the remaining accumulates in some central reservoir. The thin dotted lines are the results of MCD91 but, for 
clarity, their total rates are omitted. Our simulations were realized with 10^ super-stars. The small-scale oscillations 
present in our curves are numerical noise. See text for further comments. 
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we used 7 = 0.4 for the Coulomb logarithm. Note that 
MCD91's FP method imposes an isotropic velocity dis- 
tribution while our code allows anisotropy to develop. In 
addition to the obvious differences imposed by the use of a 
very different simulation algorithm, the following distinc- 
tions in the treatment of the physics have to be noted: 

— The collisions are treated much more realistically, on a 
particle-particle basis and outcomes are given by our 
SPH-generated grid for which realistic stellar struc- 
tures have been used. The coUisional modification of 
orbits is accounted for and mergers may occur. 

— The stell ar e volution is slightly different from MCD91 
(see Sec. 

— A "continuous" mass spectrum is used instead of the 
discrete mass classes of MCD91. To get the same av- 
erage stellar mass as these authors, the mass range is 
extended to 0.258 — 34.8 M0. Also, masses as low as 
0.01 Mq may be produced in collisions (smaller coUi- 
sional products are not allowed) while MCD91 use a 
"hard", constant minimum of 0.3 Mq. 

— We use a M*-i?, relation from MS stellar models 
(^challer et al. 1992|; [Meynct ct al. 199| ; [Charbonne] 



et al. 199£; Chabrier fc Baraffe 2000) to determine col 



MCD91 3B 



2 n 




lisional cross-sections and tidal disruption radii. 
Stellar evaporation, due to gradual energy gain 
through 2-body relaxation (see Paper I), is allowed 
in our models but MCD91 apparently enforce evolu- 
tion at constant total mass which seems reasonable be- 
cause, for a cluster with no tidal truncation, diffusive 
relaxation is expected to be inefficient. Indeed, it takes 
longer and longer to increase the (negative) energy of 
a star toward £^ > 0, as it stays for a larger and larger 
fraction of its orbital time in large-radius, low- density 
regio n s wh ere relaxation is vanishingly small ( Henon 
1960| , |1969| ). For B models with 10 super-stars, we 
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Fig. 19. Evolution of the anisotropy for 2 models of class 
B simulated with 10^ super-stars. We show the anisotropy 
parameter averaged over Lagrangian shells bracketed by 
the indicated fraction of the (remaining) stellar mass. Note 
how strong a tangential anisotropy develops in model IB, 
certainly in response to the adiabatic growth of the central 
BH. At later time, relaxation cause the central parts to 
slowly return to a more isotropic velocity distribution. The 
evolution of anisotropy in the lighter model 2B is strik- 
ingly different. For clarity, the curves have been slightly 
smoothed. 



tried to forbid relaxation-driven stellar evaporation by 
discarding "super-encounter" that lead either super- 
star to escape the system. The results appear not to 
be significantly altered by this special treatment. 

Among the results published by MCD91, those with 
which comparisons are most easily carried out and which 
are of prime interest for us, concern the growth of the 
central BH and how various processes contribute to it. In 
Fig. |lj, for all 12 models considered by MCD91, we com- 
pare the final mass of the central BH and indicate which 
process contributed most to this mass. We confirm that, 
unless the stellar mass spectrum is strongly bottom-heavy 
(case C, a = 3.5) low-mass models are dominated by stel- 
lar evolution. C models of low mass are the only ones for 
which tidal disruptions are a significant fuel source. At 
higher (initial) stellar densities, collisions dominate, with 
the densest A model as an exception. The main source of 
discrepancy between our findings and those of MCD91 is 
the more minor role of collisions in our simulations. While 
it is difficult to evaluate how MCD91's use of fixed classes 
of M, translate in their coUisional gas production rate, it 
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is certain that their n — 3 polytropes models experience 
more mass loss in off-center collisions than more realistic 
stars ( Freitag & Benz 2002b , ^ and that their mass-radius 
relation lead to an overall overestimate of collision cross- 
section. A secondary source of mismatch is our different 
prescription for stellar evolution. The temporal evolution 
of the central mass (BH + gas reservoir) for all 12 models 
is depicted in figures to p7[ 

A more detailed comparison is realized for models of 
class B for which MCD91 published the curves of the rate 
of gas production through each process. Our results are re- 
ported on_K^_^^JIer£a£aJiij_w^ 



ited number of SPH simulations by M. Davies to find fit- 
ting formulae for their outcome and incorporate collisions 
in cluster models. It is, however, doubtful that these re- 
sults, obtained with polytropic stellar models and from a 
relatively small domain of the parameter space can be ap- 
plied for realistic stars and other relative velocities and/or 
impact parameters (Freitag & Benz 2002b). 



The second important feature of the dynamics of a 
galactic nucleus, as compared to a globular cluster, is the 
likely presence of central BH with a mass in excess of 
10^ Mq (although some gl obular clusters, like M 15, may 
harbor a central BH, see Gebhardt et al. 2000 and van 



ference l with MCD91 is that their collisional rate is much der Marel 200l| ). In our code, we assume the BH stays 



higher at early times. This is probably due to the presence 
of massive stars for which their assumptions about stellar 
structure and radius should lead to the most severe over- 
estimate of collisional mass-loss. In fact, in regard of how 
different (and more detailed) our treatment of collisions 
is, it is very surprising how similar the collisional gas pro- 
duction rates are at late times. The tidal disruption rates 
are very similar at early time, with the exception of the 
least dense model (4B). At later times, our tidal gas pro- 
duction rate decreases at a steeper rate (as compared to 
MCD91) for the two densest models, while the contrary is 
true for models 3B and 4B. One possible explanation for 
the lowest late-time rates in dense models is that signif- 
icant tangential anisotropy develops in the central parts 
of these clusters, probably in response to the rapid and, 
hence, nearly adiabatic, growth of the BH, a process which 
does not significantly affect the "light" clusters. This as- 
pect is illustrated in Fig. |l^. Obviously, stars on low ec- 
centricity orbits are less likely to enter the loss cone, an 
aspect of the dynamics that MCD91 could not simulate 
with their isotropic code. On the other hand, this does 
not explain why we get a higher late time disruption rate 
for the lower density clusters. 



6. Conclusions. 

6.1. Summary 

In this second paper about our Monte Carlo code for star 
cluster simulations, we have described our inclusion of 
physical processes pertaining to the dynamics of galactic 
nuclei. 

Taking advantage of the particle-based approach of the 
MC code, collisions between MS stars are treated with a 
high level of realism. The MC sampling reproduce the rate 
of collisions between stars of various masses and the dis- 
tribution of relative velocities (and impact parameters) in 
a straightforward way. The outcome of collisions are ob- 
tained by interpolation into a comprehensive database of 
results from SPH simulations ( Freitag fc Benz 2002b| ,a|). 
This is an important improvement over previous works 
that included the role of collisions in the dynamical evo- 
lution of galactic nuclei but relied on simple-minded pre- 
scriptions for the results of collisions. In the past, only 
Ranch (1999 ) has attempted to use the outcome of a lim- 



perfectly at the center (see below) and treat its contribu- 
tion to the potential as that of a Newtonian point mass. 
The neglect of relativistic effects on stel lar orbits is p roba- 
bly a good approximation, according to Ranch (199£ ) who 
concluded that they seem to have no noticeable influence 
in his simulations. The BH grows by accreting gas released 
by the stellar system through stellar evolution, collisions 
and tidal disruptions. Whole stars may also be swallowed 
if they directly plunge through the horizon. This latter 
process completely supersedes tidal disruption for MBH 
more massive than a few 10* Mq because, then, the tidal 
disruption radius is formally inside the horizon. For the 
time being, the process of tidal disruption itself is treated 
as simply as possible, by assuming complete disruption 
of every star that enters the Roche zone around the BH. 
On the other hand, we test for super-stars entering the 
so-called "loss-cone", i.e. getting onto disruption orbits, 
in a detailed way by simulating the flne-grained diffusion 
caused by relaxation on the direction of a super-star's ve- 
locity. 

Other improvements include a simple treatment of 
stellar evolution which assumes that stars transform di- 
rectly from MS to compact remnants, in a simil ar spirit 
to what has bee n done by previous in vestigators ( Norman 
fc Scoville 1988| ; |N/Iurphy et al. 199l| ). Also, we have im- 
plemented "particle doubling" to maintain high resolution 
even in simulations where a lot of stars are either destroyed 
or ejected from the cluster. 

These new features have been extensively checked 
against (semi-) analytical predictions and simulations from 
the literature. In most cases, the tests are highly success- 
ful. In particular, collision rates are nicely reproduced, not 
only when integrated over the whole cluster but also as a 
function of distance from the center and of the masses of 
stars. The effects on the stellar cluster of an adiabatically 
growing central black hole are nearly perfectly in agree- 
ment with theoretical predictions. The standard "Bahcall 
and Wolf" density cusp is obtained in the case 

tidal disruptions are taken into account but collisions are 
switched off or inefficient. In highly collisional models, a 
shallower cusp, with exponent around —0.5 is produced, in 
good agreement with what was reported in previous stud- 
ies. Gas production by the stellar cluster through various 
processes (tidal disruptions, collisions, stellar evolutions) 
are also in good agreement with results from the litera- 



M. Freitag & W. Benz: Monte Carlo Cluster Simulations II. 



27 



ture, obtained with a variety of numerical methods. Most 
of the discrepancies can be easily explained. In particu- 
lar, it appears that the role of collisions has been over- 
estimated in previous works, due to over-simplified as- 
sumptions about the coUisional outcome (complete disrup- 
tions or simple semi-analytical treatment applied to poly- 
tropic models) and, maybe, to their being included into 
the simulations in a quite nonphysical way, in the case 
of direct Fokker-Planck methods. Concerning tidal dis- 
ruptions, some disagreement, for which we have found no 
straightforward explanati on, is observed with the works of 
Duncan & Shapiro (1982 ) and Amaro-Seoane fc Spurzcm 



(200l|)J These mismatches are not se vere, however, and 



as the resolution of the simulations by Duncan fc Shapiro 
(19821)" was quite lowp| and the results plotted by Amaro 
Seoanc fc Spurzem (2001 ) come only from preliminary 



computations, we can not draw definitive conclusions from 
these comparisons. Furthermore, there is no clear trend in 
these differences and we get better agreements in other 
cases (with, e.g., model I of Duncan fc Shapiro 1983| ), a 
fact which seems to exclude any important flaw in our al- 
gorithm. Unfortunately, iV-body methods seem still a long 
way from allowing simulations of the relaxational dynam- 
ics around a black hole and, thus, providing more direct 
check of our approach and, more generally, of the appli- 
cability of the loss-cone theory and the Chandrasekhar 



treatment of relaxation in such a situation (e.g., Spurzem 

fc Kug fsl 2ooq ). 



6.2. Future work 

In Sec. 8.2 of paper I, we have already mentioned many 
improvements/additions that we plan to incorporate in 
future versions of the code. Here, we update and complete 
this list: 



1. 



Capture of compact stars by the central BH 
through emission of gravitational radiation. This 



process has been presented in Sec. |l.l| . Predicting the 
rate and characteristics of these events has recently 
become a main focus of our work and very encouraging 



results have already been reported in Freitag (2001). 
Refined treatment of stellar evolution. The most 
severe shortcoming of our present modeling of SE is the 
absence of giant phase. Taking it into account should 
greatly enhance the number rates of collisions and tidal 



stripping ( 


Davies et al. 1998; 


Bailey & Davies 1999; 


Magorrian fc Trcmaine 1999; 


Syer fc Ulmer 1999|) al- 



though the amount of released gas may be limited due 
to the very low density of giants' envelopes and this 
may not increase the BH's growth as this gas would 
be liberated anyway through stellar evolution. Others 
aspects of SE that we shall incorporate are: progres- 
sive mass loss on the MS, natal kicks for neutron stars 
and coUisional rejuvenation. 



They used a few thousands particles but their cloning al- 
gorithm increased the relative resolution at large negative en- 
ergies, i.e. close to the BH. 



Refined treatment of tidal interactions. 

We should treat the hydrodynamical nature tidal 
disruptions with the same level of realism that we 
achieved for collisions. This will be essential if we 
want to cope with envelope-stripping of giant stars ( pi| 
Stefano et al. 2001) and other "tidally perturbed" stars 
(Alexander & Livio 2001). Stars can also be tidally 
captured by the central BH. As more and more or- 
bital energy is transfered to oscillations at each subse- 
quent pericenter passage, disruption is the most prob- 



able outcome (Novikov et al. 1992). 
Assuming that the BH stays fixed at the center of the 
cluster is an over-simplification. If the central BH's 
wandering is of larger extent than its tidal disruption 
radius i?disr, there will be no regime of empty loss cone 
( Bigurdsson fc Rees 1997 ) . For a cluster with core ra- 
dius Rc, equipartition predicts a wandering radius of 
order » Rc\/M*/Mbk (|BahcaU fc Wolf 1976| ; ^ 
fc Tremaine 1980| ; [Chatterjee et al. 2002] ), and 



Ry. 
Rdi! 



400 



Rc 
1 pc 



i?0 



Mo 



106 Mp 



See Magorrian fc Trcmaine (1999| ) for hints at the 



possible effects of the wandering on the tidal disrup- 
tion rate. Young (1977) made a rough estimate of the 
correction and deemed it not to alter the disruption 
rate drastically. However, as suggested by Alexander 
fc Livio (2001[ ), these motions of the BH may allow 
stars that have been tidally perturbed to escape fur- 
ther, disruptive, close interactions with the BH, which 
is of high potential interest for the Galactic center. 
Large angle scatterings. 

2-body gravitational encounters with impact parame- 
ter of order or smaller than bo — G{Mi + M2)/ V^^j lead 
to scattering angles of order it. Although they only con- 
tribute a fr action ln(5m ax/&o)~^ < 0.1 to the overall 
relaxation ( Hcnon 1973| , p. 198), the y may dominate 
the rate of evaporation from the cusp (Lin fc Tremaine 
198C ; [Goodman 1983|) and of c aptures on relativistic 
orbits ( Bigurdsson fc Rees 1997 ) . Such "kicks" can not 



be decomposed into smaller defiections but can proba- 
bly be introduced explicitly in the MC code in a similar 
way as collisions. 
Inclusion of binary stars. 



In a "normal" population (Duqucnnoy fc Mayor 1991 



most binaries have (internal) orbital velocities smaller 
than the velocity dispersion near the central BH (a 
few hundreds kms~^) and will eventually be disrupted 
through interactions with other stars. However, some 
small fraction may be hard enough to survive and 
evolve into compact binaries. Whether hard binaries 
will have an important dynamical role has to be ex- 
plored. Their interaction with the central BH is of par- 
ticular interest. Indeed, if it passes sufficiently close to 
the BH, a binary will be tidally disrupted with the 
likely result of ejecting one star out of the cluster at 
very high velocity and leaving the other one bound to 
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the BH pills 1988| , |1991| ). This is another channel to 
form extreme mass ratio binaries to be detected by 
LISA. 

Interaction with a central accretion disk or gas 
cloud. 

The early evolution of galactic nuclei may well lead 
to the accumulation of a quasi spherical central gas 
cloud with high enough a density to interact strongly 
with the stellar cluster. This situation has not yet been 



given the attention it deserves (see, however, Bpurzem 
1992, and references therein) but further investiga- 



laborators ( 


Amaro-Seoane & Spurzem 2001 


Seoane et al. 2001 


)• 



In AGNs, stars may be captured by an accretion 
disk through repeated impacts which can strongly re- 
shape the stellar distribution in the vicinity of the 



BH (iNorman & Silk 198S; |Syer et al. 1991; 


Ranch 


1995|; Vokrouhlicky & Karas 199^; Karas & Subr 2001; 


Vilkoviskij & Czerny 2002 


). The further stellar and 



orbital evolution of the disk-embedded stars is a com- 
plex subject. Interesting possibilities include enhanced 
rate of collisions and growth of massive stars by accre- 
tion of disk material. Note that even if the interactions 
with the accretion disk are not efhcient enough to grind 
down orbits into the disk, stellar formation probably 
occurs in situ (Goodman 2002) so that the presence of 



stars in the disk has to be expected anyway. A possi- 
ble way of accounting for the role of the accretion disk 
in numerical models would be to use the MC code to 
simulate the outer quasi-spherical parts of the clus- 
ter where relaxati on is import ant and couple it with 



a code like that of Subr (2001) which treats the inner 



regions, where interactions with the disk dominate the 
dynamics, in axisymmetrical geometry. 
Gas dynamics. 

Including stellar evolution without a better prescrip- 
tion for the fraction of gas that eventually finds its way 
to the central BH is nearly pointless, as demonstrated 



by simulations in 


Freitag (200C). Early studies ( 


Bailey 


1980; 


Loose feFricke 1980; 


David et al. 1987a 


b; 


Kunze 


et al. 1987 


; Norman & Scoville 1988) concluded that 



most of the gas finds its way to the central BH but 
they lacked detailed account of the feed-back on the 
gas of the energy released by the central source and su- 
pernova explosions and of the complex, non-spherical, 



evolving geometry of the gas flow (see, e.g., Williams 



et al. 1999; Ciotti & Ostriker 2001, for recent attempts 
at tackling these intricacies). 

This list can be lengthened virtually without end. But 
before we hurry and include more and more complexity in 
our simulations, we must keep in mind that each new pro- 
cess to be added comes with its own uncertainties of both 
physical and numerical nature, so that the impression of 
added "realism" may be misleading. In such a context, it 
is all the more useful to dispose of a numerical tool flex- 
ible enough to allow changes in the treatment of various 



physical effects and fast enough to allow large sets of sim- 
ulations to be conducted to test for the influence of these 
modifications. 

Another line along which we have to progress is to 
develop definite observational predictions. Here are a few 
examples: 

— Surface luminosity and color profiles for central cusps. 

— Rate and characteristics of radiation fiarcs following 
the tidal disruption of a star. 

— Appearance (and radial distribution) of stars modified 
by collisions or tidal interactions with the MBH. 

— Rate and characteristics of gravitational waves signals 
from captured stars. 

All examples but the first are complex problems of 
their own and have already been the subject of many de- 
tailed, if not conclusive, studies. Fortunately these aspects 
are essentially decoupled from the cluster dynamics, in the 
sense that they have no obvious back-influence on it, so 
that we should be able to "map" results from the literature 
on the outcome of our simulations. 
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Appendix A: Building of initial models of galactic 
nuclei 

To obtain initial cluster realizations for our simulations, 
we proceed in two stages: (1) We set the radii Ri, specific 
kinetic energies, Ti and moduli of specific angular mo- 
mentum, J,j of all super-stars^ while trying to ensure dy- 
namical equilibrium. (2) We set the stellar masses of the 
super-stars, M*, according to a given initial mass func- 
tion (IMF). To get an aged stellar population, we may 
also evolve this IMF according to the "ZAMS — >remnant" 
relation specified in Sec. 4.1. As the number of stars a 
super-star stands for must be the same for all super-stars, 
this stage also implicitly determines the super-star's mass. 
Mi = {Nf,/Np)M* where A'p is the number of super-stars 
the model consists of and A^, is the number of stars rep- 
resented by the model. 



http : //obswww.unige . ch/~pf ennige/gravitor/ 
gravitor_e .html 

Remember that a super-star actually represents a spherical 
shell of stars. 
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A.l. Positions and velocities 

The safest way to obtain a system that is not only virial- 
ized {2Tc\ + Uc\ = where Td is the total kinetic energy 
and Uci the total gravitational energy), but a genuine sta- 
tionary solution of the coUision-lcss Boltzmann equation, 
is to start from a onc-particle DF /(X, V) which depends 
on the position X and velocity V only through isolating 
integrals of motions, namely E and J, for a stellar cluster 
that obeys spherical symmetry (Binney & Tremaine 1987, 
Chap. 4), 



(2) For each super-star, radius Ri is randomly selected ac- 
cording to the probability density dMr/dR. This is 
done by creating a random number Xj-an with uniform 
probability over [0; 1[ and (numerically) inverting the 
MriR) relation: Ri = M-i(Xi.a„). 

(3) Once the radius Ri of super-star i is determined, we 
have to select a velocity Vi according to distribution 
g{V) (X V^Fj^V ^ + <^{Ri)). Here we use a simple re- 

Sec. 



fiX,V)^FiEiX,V),J{X,V)), 



(A.l) 



with 



EiX,V) = -V^ 



$(i?) and J{X,V) = RV±, (A.2) 



where R= \X\, V — \V\, V± is the modulus of the compo- 
nent of V perpendicular to X (with the cluster center as 
origin of coordinates) and <I> is the (smooth) gravitational 
potential. For the sake of simplicity, we only considered 
initial cluster models with isotropic velocity distributions 
for which F is a function of E only. Note that our MC 
code can tackle any velocity distribution and that some 
level of anisotropy develops during the run of most cluster 
simulations. 

$ is itself determined by the DF through Poisson equa- 
tion: 



d$ 



R- 



d2$ 



'dR ' ' di?2 
with the density p given by 

p(<i>) = in 



47rG/9($), 



dVV^F(-V^ 
^2 



$). 



(A.3) 



(A.4) 



It is customary to define so-called relative energy and 
potential through 



-r dcf 



$ and 



dcf 

e = $0 



E 



(A.5) 



with $0 chosen so that F{£) = for e < 0. For a cluster 
of finite radius i?ci, — —GMd/Rd- 

Thus, to build a cluster model, we do the following: 

(0) Choose an expression for F{e). Traditional choices are, 



among others, Plummer's or King's models (Binney & 



Tremaine 1987). 



(1) Integrate *(i?) and Mr{R) with a Runge-Kutta 
scheme ( Hairer et al. 198"7| ): 

*d = -47rGp(*) - 2 ' 
Mr J \ 47rp(*)i?2 




(A.6) 



jection method (Press et al. 1992 



7.3) with a 

constant upper bound given by — 2$(i?j)F(<i>(i?i)).[^ 
The specific kinetic energy of the super-star is thus 
Ti = To set the specific angular momentum Ji 

with account of isotropy, we generate another random 
number Xian and compute Vmd — ^^i(l — 2Xran) and 

Jr=RWV'-yrL- 

(4) Finally, perfect virial energy balance is enforced by a 
slight re-scaling of the velocities. 



Each evaluation of the function p(^') requires itself a 
numerical integration of Eq. A.4. The integration of 
system |A.6| is terminated either when the relative po- 
tential reaches (for tidally truncated models) or when 
Mr has attained some asymptotic value. At that point, 
we have obtained array representations of i?, 5", p and 
Mr- We re-normalize them to the "TV-body" system of 
units (see Sec. O). 



In its present form, this procedure does not explicitly 
allow for a central BH. But if we add such a point mass at 
the center with a very small mass (as compared to Md), 
it will only slightly perturb the potential energies of the 
innermost super-stars and the resulting system will still 
be very close to dynamical equilibrium. This is the reason 
why we must always start simulations with "seed" black 
holes instead of already grown (super-)massive ones. An 
advantage of this method is that the integrated influence 
of the BH's growth on the stellar system is "automati- 
cally" computed! The main drawback is that we cannot 
start with models that represent today's galactic nuclei 
but have to guess initial conditions that lead to such con- 
figurations after a Hubble time. This has not yet been 
explored systematically. 

The cluster produced with this algorithm has no mass 
spectrum, i.e. all super-stars have the same mass Mp — 
Mc\/Np. We now explain how we construct a stellar mass 
spectrum. 

A.2. asses 

We model IMFs that are piece- wise power-laws, 
dA^. 



dM, 



oc m: 



for Mk-i < M^ < Mk, 



(A.7) 



between some Mq — M,„in and AIk — -^max- 

For a given set of Mj, (k — 0, . . . , K) and 
(fc = 1, . . . , K). The un-normalized number of stars with 
masses < Af, is, for Mk-i < Af* < Mk'. 



NiM,) = Nk-i+Ck I TT^dM* 



dM, 



1 



k-l 



(M, 



k-l 



(A.8) 



1 - Ctk 

with Ck = Cfc-iM^"""""-'^ (we can set Ci = 1). Once 
the Nk have been computed, we randomly determine the 

Bound particles have /2 + ^{R) < 0. Furthermore, well- 
behaved DF have dF/dE < so that the maximum value at a 
given R is F{${R)). 
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stellar mass of each super-star in turn. We first generate 
a random number iVran with uniform [0; Nk] distribution 
{Nk is the un-normalized total number). We then find 
index j such that iVj_i < A^ran < Nj and invert N{M^) 
to find the stellar mass for super-star i: 

M* = Im}:^' + (1 - aj) ) ~ ■ (A.9) 

Note that we never need to state the actual total num- 
ber of stars (or, equivalently, the total mass in Mq) or 
the size of the cluster in pc when building initial models. 
This must only be specified before starting an evolution- 
ary Monte Carlo simulation as these mass and size scales 
determine the relative importances of various processes 
(e.g. relaxation vs. collisions) and allows to translate the 
A''-body time units into years. 
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